MateWiki - Contribuciones del usuario [es]

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-10T09:54:14Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

{{ TrabajoED | Flujo alrededor de un obstáculo circular. Grupo 68| [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] | Nicolás Rus Pérez, Jeremy García Herrera, Antonio Rodriguez Montes, Franco Luis Santillan, Juan Jose Salhua Palma.}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]
Un fluido es un medio contínuo deformable con la propiedad de que no transmite tensiones tangenciales. Estos no tienen una forma definida y por tanto adoptan la forma del recipiente en el que se encuentran, y pueden encontrarse tanto en estado líquido (agua) como gaseoso (aire).

En este trabajo se estudiará el flujo de un fluido alrededor de un obstáculo circular considerando a este un fluido incompresible. Esta propiedad indica que el fluido no disminuye su volumen al serle aplicado una fuerza, y es propia de los fluidos en estado líquido.

Por conveniencia se trabajará en coordenadas con base cilíndrica (polares).

==Mallado==
Se genera a continuación un mallado en el que se representan los puntos interiores a la región ocupada por el fluido. Por convenio se representará el obstáculo circular centrado en el origen de coordenadas y con radio 1. De esta forma la región ocupada por el fluido será todo el exterior de la circunferencia unidad.

El objetivo de realizar este mallado es dividir el medio continuo (la región ocupada por el fluido) en subdominios más pequeños denominados elementos finitos. El conjunto de estos elementos forman el mallado que cubre toda la región y permite realizar cálculos y simulaciones numéricas particulares en cada uno de ellos.

Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo <math>[-4,4]×[-4,4]</math>.
[[Archivo:Mallado1g68.JPG|thumb|right|Región del fluido y obstáculo]]
{{matlab|codigo=

r=linspace(1,5,50); %Definen el rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Genera el mallado

hold on
X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);
Z=0.*R;
mesh(X,Y,0*Z); %Dibuja la región
plot(cos(a),sin(a),'k','lineWidth',1); %Dibuja el obstáculo
axis([-4 4,-4 4]); %Se fijan los ejes
title ('Regíon Fluido'); %Definen los ejes y el nombre de la gráfica
xlabel 'EJE X'
ylabel 'EJE Y'
hold off
}}

==Campo de velocidad de las partículas==
Una vez definido el mallado de trabajo del fluido y el obstáculo, se procede a calcular las velocidades de las partículas del fluido en los distintos puntos de la región ocupada por el fluido incompresible. Partimos sabiendo que el campo de velocidades viene definido por el gradiente de la siguiente función potencial:

<center><math>\varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) </math></center>

A continuación se representa graficamente la función potencial:
[[Archivo:Mallado222g68.JPG|thumb|400px|right|Función potencial en la región del fluido]]
{{matlab|codigo=

r=linspace(1,5,30); %Define rho (r) y theta (a)
a=linspace(0,2*pi,30);

[R,A]=meshgrid(r,a); %CGenra el mallado

hold on
X=R.*cos(A); %Parametrización la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Define la función potencial
Z=f(R,A);
surf(X,Y,Z); %Representa la función
plot(1*cos(a),1*sin(a),'k','lineWidth',3); %Representa el obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de colores, define el titulo de la gráfica y los ejes
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
hold off
}}
Para calcular el campo de velocidades de las partículas del fluido es necesario el campo gradiente de la funcion potencial anterior, y dicho gradiente viene definido por la siguiente expresión:

<center>\(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)</center>

Aplicando dicha expresión a la función potencial, el gradiente resultante sera:

<center><math> \vec u(\vec e_\rho,\vec e_\theta,\vec e_z)=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta)]\vec e_\theta </math></center>

Dado que la representación en MATLAB se realiza en coordenadas cartesianas <math>(\vec i,\vec j,\vec k)</math>, y el gradiente calculado se encuentra en coordenadas ciliíndricas <math>(\vec e_\rho,\vec e_\theta,\vec e_z)</math> ,debemos de pasar el gradiente de la función potencial a las primeras. Esto es posible gracias a la matriz de cambio de base entre coordenadas cilíndricas y cartesianas, por lo que la multiplicamos por el vector gradiente dado <math>(\vec u)</math> en coordenadas cilíndricas.

Matriz de cambio de base de coordenadas cilindricas a coordenadas cartesianas:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z\\
\end{pmatrix}
</center>

Por lo que:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
(1-\frac{1}{\rho^2})\cos(\theta)\\
-(1+\frac{1}{\rho^2})\sin(\theta)\\
0\\
\end{pmatrix}
</center>
Resolviendo a coordenadas cartesianas:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta)\right) \vec j </math></center>

A continuación representamos el campo de velocidades de las partículas sobre las superficies equipotenciales del flujo del fluido.
[[Archivo:Mallado3g68.JPG|thumb|500px|right|Campo de velocidades de las partículas sobre las líneas equipotenciales del fluido]]
{{matlab|codigo=

r=linspace(1,5,50); %Define rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Crea la malla

X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Definimos la función potencial
Z=f(R,A); %Aplica la función
contour(X,Y,Z,30); %Dibuja las curvas de nivel
hold on %Define las componentes X e Y del gradiente
Cx=(1-(1./R.^2)).*cos(A).^2+(1+(1./R.^2)).*sin(A).^2;
Cy=(1-(1./R.^2)).*sin(A).*cos(A)-(1+(1./R.^2)).*sin(A).*cos(A);
quiver(X,Y,Cx,Cy); %Dibuja el campo de velocidades
plot(1*cos(a),1*sin(a),'k','lineWidth',1); %Representa del obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de color, titulo y nombre a los ejes
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}
Gracias a la gráfica anterior se observa como los vectores de velocidad definidos por <math> \vec u </math> son siempre perpendiculares a las líneas equipotenciales del fluido.
<center>[[Archivo:Mallado33.JPG]]</center>

==Rotacional Nulo y Divergencia Nula==
===Rotacional Nulo===
El rotacional nulo del campo de velocidades indica que el flujo es '''irrotacional''' , es decir , carece de vórtices o giros locales. Esto significa que una partícula de fluido no rota sobre sí misma mientras se desplaza, característica típica de los flujos potenciales.

El operador rotacional nos indica la dirección y la velocidad del giro del campo en cada punto.

A continuación se demuestra que el rotacional es nulo:

<center><math>
\nabla\times \vec{u} = \frac{1}{\rho} \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & \rho(-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix}=</math>

<math>
= \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} \left[ \rho \left( -(1 + \frac{1}{{\rho^2}})\sin(\theta) \right) \right]\vec{e}_z - \frac{\partial}{\partial \theta} \left[ \cos(\theta) \left(1 - \frac{1}{\rho^2}\right) \right] \vec{e}_z \right]=\frac{1}{\rho} \left[ -(1 - \frac{1}{\rho^2})\sin(\theta) + \sin(\theta)(1 - \frac{1}{\rho^2}) \right]=0
</math>

<math>\boxed { \nabla \times \bar { u } =0 } </math></center>

Se obtiene que el rotacional es 0, lo cual indica que las partículas del fluido no están girando

===Divergencia Nula===
Interpretación de la divergencia nula: La condición ∇·u = 0 corresponde a un flujo incompresible, donde la densidad del fluido permanece constante y el volumen de cualquier elemento material no cambia a lo largo del tiempo. Esto refleja la conservación local de la masa en ausencia de fuentes o sumideros.

<center><math>
\nabla \cdot \vec{u}(\rho, \theta, z) = \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} (\rho u_\rho) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (\rho u_z) \right\}= </math>

<math>
= \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} \left( (\rho - \frac{1}{\rho})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) \right) \right\}=
</math>
<math>
\frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0
</math>

<math>\boxed { div(\bar { u } )=0 } </math></center>

== Líneas de corriente ==
Las líneas de corriente del flujo alrededor de un cilindro circular (en sección transversal) presentan un patrón simétrico con respecto al eje horizontal (eje x) y antisimétrico con respecto al eje vertical (eje y).
[[Archivo:Linea_corrienteTangencial.png|thumb|400px|right|Líneas de corrientes]]

Las líneas de corriente <math>\psi = \text{cte}</math> describen las trayectorias del flujo potencial alrededor del obstáculo circular, mostrando simetría y adaptación suave a la geometría. Estas líneas son tangentes al campo de velocidades en cada punto, verificando que <math>\vec{u}\cdot\nabla\psi = 0</math> . El patrón obtenido revela aceleración del fluido cerca del cilindro, condición de no penetración en la superficie, y recuperación del flujo uniforme lejos del obstáculo, características típicas de un flujo incompresible e irrotacional.

<math>\vec{v}</math> se calculará con : <math>\vec{v}=\vec{k}\times \vec{u}</math>

<math>\vec{v}</math> es irrotacional dadon que el rotacional calculado en el apartado 3.1 es nulo.
De igual forma iene un potencial escalar <math>\psi </math> el cual tiene <math>\vec{v}</math> de gradiente.

Para el potencial escalar <math>\psi </math> se resolverá este sistema de ecuaciones: <math>\nabla\psi =\vec{v}</math> cumpliendo la condición anterior de líneas <math>\psi =cte</math>.

Empezamos las operaciones:
<math>\vec{v} =\vec{k}\times \vec{u} = \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
0 & 0 & 1 \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & (-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

<center><math>\boxed{\frac{\partial \psi}{\partial \rho} = V\rho,
\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = V\theta,
\frac{\partial \psi}{\partial z} = Vz}
</math></center>

<math>\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho = \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

<math>\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})\Rightarrow\psi(\rho,\theta,Z)= \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

Dando como resultado
<math>\psi: \sin{\theta}(\rho-\frac{1}{\rho})=cte</math>
[[Archivo:Psi_lineasCorriente.png|thumb|400px|right|Líneas de corrientes]]
{{matlab|codigo=
% Definir dominio en coordenadas cilíndricas
r = linspace(1, 5, 50); % Coordenada radial (ρ) desde r=1 (obstáculo)
a = linspace(0, 2*pi, 50); % Coordenada angular (θ)

% Crear malla
[R, A] = meshgrid(r, a);

% Convertir a coordenadas cartesianas
X = R .* cos(A);
Y = R .* sin(A);

% 1. FUNCIÓN POTENCIAL φ (phi)
phi = (R + (1./R)) .* cos(A);
% 2. CAMPO DE VELOCIDADES u = ∇φ

% Componentes en coordenadas cilíndricas:
% u_ρ = (1 - 1/ρ²)cos(θ)
% u_θ = -(1 + 1/ρ²)sin(θ)

u_rho = (1 - (1./R.^2)) .* cos(A);
u_theta = -(1 + (1./R.^2)) .* sin(A);

% Convertir a coordenadas cartesianas
u_x = u_rho .* cos(A) - u_theta .* sin(A);
u_y = u_rho .* sin(A) + u_theta .* cos(A);

% 3. FUNCIÓN DE CORRIENTE ψ (psi)
psi = (R - 1./R) .* sin(A);

% 4. GRÁFICAS

figure('Position', [100, 100, 1200, 800]);

% Graf 1: Líneas de corriente (ψ) y campo de velocidades
subplot(2, 2, 2);
contourf(X, Y, psi, 30, 'LineStyle', '-', 'LineWidth', 0.5);
hold on;
colorbar;
title('Líneas de Corriente (ψ) y Campo de Velocidades');
xlabel('EJE X');
ylabel('EJE Y');
axis equal;
axis([-4, 4, -4, 4]);
grid on;

% Dibujar campo de velocidades
quiver(X(1:3:end, 1:3:end), Y(1:3:end, 1:3:end), ...
u_x(1:3:end, 1:3:end), u_y(1:3:end, 1:3:end), ...
1.5, 'w', 'LineWidth', 0.7, 'MaxHeadSize', 0.5);

% Dibujar obstáculo circular
plot(cos(theta_circ), sin(theta_circ), 'k', 'LineWidth', 2, ...
'Color', [0.2, 0.2, 0.2]);
fill(cos(theta_circ), sin(theta_circ), [0.8, 0.8, 0.8], ...
'EdgeColor', [0.2, 0.2, 0.2], 'LineWidth', 1.5);
hold off;
}}

== Velocidad extrema y puntos de remanso==

En la frontera del obstáculo, definida por la circunferencia unidad <math>\rho=1</math>, el módulo de la velocidad viene dado por

<center><math>|\vec u| = 2|\sin\theta|</math></center>

De esta expresión se deduce que la velocidad máxima se alcanza cuando

<math>\sin\theta=\pm1</math>, es decir, en

<center><math>\theta=\frac{\pi}{2},\ \frac{3\pi}{2}</math></center>

que corresponden a los puntos cartesianos

<center><math>(0,1)\ \text{y}\ (0,-1)</math></center>

Estos son los puntos donde la velocidad del fluido es máxima sobre la frontera.

Por otro lado, la velocidad se anula cuando

<center><math>\sin\theta=0 \Rightarrow \theta=0,\ \pi</math></center>

lo que corresponde a los puntos

<center><math>(1,0)\ \text{y}\ (-1,0)</math></center>

Estos puntos reciben el nombre de '''puntos de remanso'''.

Para verificarlo numéricamente en MATLAB se ha representado el módulo de la velocidad sobre la frontera:
[[Archivo:Subir1.png|thumb|500px|center|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular. La velocidad es nula en los puntos de remanso y máxima en la dirección vertical.]]
{{matlab|codigo=
nr=100; nt=200;
r=linspace(1,5,nr);
th=linspace(0,2*pi,nt);
[R,TH]=meshgrid(r,th);

ur=(1-1./R.^2).*cos(TH);
uth=-(1+1./R.^2).*sin(TH);

Ux=ur.*cos(TH)-uth.*sin(TH);
Uy=ur.*sin(TH)+uth.*cos(TH);

Umag=sqrt(Ux.^2+Uy.^2);

UmagS=Umag(:,1);
THS=TH(:,1);

[Umax,imax]=max(UmagS);
[Umin,imin]=min(UmagS);

thetaMax=THS(imax);
thetaMin=THS(imin);

figure
plot(THS,UmagS), grid on
hold on
plot(thetaMax,Umax,'ro',thetaMin,Umin,'bo')
title('Modulo de la velocidad sobre la frontera S')
xlabel('\theta'), ylabel('|\itu|')
hold off
}}

== Presión del fluido y ecuación de Bernoulli==

Se supone que la densidad del fluido es constante e igual a

<center><math>\rho=2</math></center>

y se cumple la ecuación de Bernoulli

<center><math>\frac12\rho|\vec u|^2+p=C</math></center>

Despejando la presión se obtiene

<center><math>p=C-\frac12\rho|\vec u|^2</math></center>

Tomando para la constante el valor <math>C=1</math>, la presión queda completamente determinada en toda la región del fluido.

El campo de presiones se ha representado en MATLAB mediante:

{{matlab|codigo=
rho=2;
C=1;

P=C-0.5*rho*Umag.^2;

PS=P(:,1);
[Pmax,ipmax]=max(PS);
[Pmin,ipmin]=min(PS);

thetaPmax=THS(ipmax);
thetaPmin=THS(ipmin);

figure
contourf(R.*cos(TH),R.*sin(TH),P,30,'LineStyle','none');
colorbar
hold on
viscircles([0 0],1,'Color','k')
axis equal
title('Presión p(x,y)')
xlabel('x'), ylabel('y')
hold off
}}
{| class="wikitable" style="border:0; margin:auto; background:none;"
| style="border:0;" | [[Archivo:Subir3.png|thumb|500px|Campo de presiones p(x,y) según la ecuación de Bernoulli]]
| style="border:0;" | [[Archivo:Subir2.png|thumb|500px|Módulo de la velocidad ∣𝑢⃗∣ sobre la frontera del obstáculo 𝑟=1 en función del ángulo 𝜃]]
|}
Se observa que la presión es máxima en los puntos de remanso, donde la velocidad es nula, y mínima en los puntos donde la velocidad alcanza su valor máximo. Esto concuerda con la ecuación de Bernoulli, ya que un aumento del módulo de la velocidad implica necesariamente una disminución de la presión.

==Velocidad y presión al rodear el obstáculo.==

Ya teniendo el campo de presiones se puede apreciar que tiene relación con la velocidad, ya que cuando disminuye la presión aumenta el módulo de la velocidad, tratándose de un resultado lógico debido a que cuando hablamos de un fluido incompresible su valor de Trinomio de Bernoulli ha de ser constante.
Lo podemos ver mediante la siguiente expresión: (Ec. Bernoulli).

<math>B=z+P/γ+v^2/2g</math>

B=valor del trinomio de Bernoulli

P= presión

z= altura

v=velocidad que lleva el fluido en determinado punto

g=valor de la gravedad

γ=valor del peso específico

Donde podemos ver que si la presión aumenta la velocidad tiene que disminuir y viceversa.

==Paradoja de D'Alembert.==

La magnitud pn representa la fuerza que el fluido ejerce de forma puntual sobre la superficie del obstáculo. Para estudiar si existe empuje neto en la dirección horizontal, basta analizar la suma de las proyecciones de todas esas fuerzas sobre el vector i. Si dicha suma resulta nula, concluiremos que no hay fuerza total en esa dirección.

Del apartado anterior conocemos ya la presión en la frontera del círculo unidad:
<math>p=10-∥u∥^2=9-4 sin⁡^2 θ-4sin⁡θ</math>

Cálculo del vector normal

Sobre la circunferencia de radio 1 el vector normal exterior es

<math>n=-eρ.</math>

Para expresarlo en coordenadas cartesianas utilizamos la matriz de cambio de base de polares a cartesianas. Tomando su traspuesta se obtiene:

\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}

Aplicando esta matriz a <math>(1,0,0)^T</math>, que corresponde al versor radial eρ, se obtiene:
\begin{pmatrix}
\cos \theta &\\
sen \theta &\\
0
\end{pmatrix}
<math>=cos⁡θi+sin⁡θj</math>

Al tomar el signo negativo para el normal exterior,

<math>n=cos⁡θi-sin⁡θj</math>,

de donde su proyección sobre i es:

<math>n⋅i=-cos⁡θ</math>.

Integración sobre la frontera

El borde del dominio se parametriza como
<math>(ρ,θ)=(1,θ), con θ∈[0,2π)</math>.

La contribución total de la presión en la dirección horizontal es:
<math>\int_{0}^{2π}p(n⋅i)dθ=\int_{0}^{2π}p(-cos⁡θ)dθ.</math>

Sustituyendo la expresión de p:
<math>\int_{0}^{2π}(9-4sin⁡^2 θ-4sin⁡θ)(-cos⁡θ)dθ.</math>

Integramos término a término:
<math>\int_{0}^{2π}-9cos⁡θdθ+\int_{0}^{2π}4sin⁡^2 θcos⁡θdθ+\int_{0}^{2π}4sin⁡θcos⁡θdθ.</math>

Cada una de estas integrales es periódica con una primitiva continua en todo el intervalo; evaluándolas entre 0 y 2π se obtiene:
<math>[-9sin⁡θ]_{0}^{2π}=0,[(4sin⁡^3 θ)/3]_{0}^{2π}=0,[2sin⁡^2 θ]_{0}^{2π}=0.</math>

Por tanto,
<math>\int_{0}^{2π}p(n⋅i)dθ=0.</math>

Conclusión

La componente horizontal de la fuerza total ejercida por el fluido sobre el obstáculo es nula. En consecuencia, el fluido no produce empuje neto en la dirección i.

==Cambio de función potencial==

En este apartado repetimos los procesos de los puntos (2), (3) y (4), pero tomando ahora como función potencial

<center><math>\varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}</math></center>

El obstáculo sigue siendo la circunferencia unidad centrada en el origen, y el fluido ocupa la región exterior, trabajándose de nuevo en coordenadas cilíndricas (polares).

===Representación de la función potencial y del campo de velocidades===

En primer lugar representamos la nueva función potencial <math>\varphi</math> sobre la región ocupada por el fluido, manteniendo el mismo mallado del apartado inicial. Esta representación permite visualizar cómo el término angular <math>\theta/(4\pi)</math> introduce una ligera inclinación en la superficie potencial respecto al caso anterior.

[[Archivo:nuevafuncionpotencial.png|thumb|400px|right|Nueva función potencial del flujo con término angular]]

La velocidad de las partículas del fluido viene dada de nuevo por el gradiente del potencial. En coordenadas cilíndricas, el campo de velocidades es

<center><math>\vec u(\rho,\theta,z)=\nabla\varphi
= u_\rho\,\vec e_\rho+u_\theta\,\vec e_\theta+u_z\,\vec e_z</math></center>

donde

<center><math>
u_\rho=\frac{\partial\varphi}{\partial\rho}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta, \qquad
u_\theta=\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}
= -\left(1+\frac{1}{\rho^2}\right)\sin\theta+\frac{1}{4\pi\rho},\qquad
u_z=0.
</math></center>

Para poder representarlo en MATLAB se pasa el campo a coordenadas cartesianas mediante la matriz de cambio de base de cilíndricas a cartesianas:

<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k
\end{pmatrix}
=
\begin{pmatrix}
\cos\theta & -\sin\theta & 0\\
\sin\theta & \cos\theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z
\end{pmatrix}
</center>

Aplicando esta matriz al vector <math>(u_\rho,u_\theta,0)^T</math> se obtiene

<center><math>
\vec u(\vec i,\vec j,\vec k)
=\left[(1-\tfrac{1}{\rho^2})\cos^2\theta+(1+\tfrac{1}{\rho^2})\sin^2\theta
-\frac{\sin\theta}{4\pi\rho}\right]\vec i
+\left[(1-\tfrac{1}{\rho^2})\sin\theta\cos\theta
-(1+\tfrac{1}{\rho^2})\sin\theta\cos\theta
+\frac{\cos\theta}{4\pi\rho}\right]\vec j.
</math></center>

Sobre las curvas de nivel de <math>\varphi</math> se representa este nuevo campo de velocidades:

[[Archivo:campodevelocidadesnuevopotencial.png|thumb|500px|right|Campo de velocidades sobre las líneas equipotenciales para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.1: Potencial y campo de velocidades para la nueva función φ

r = linspace(1,5,30); % Coordenada radial (rho)
a = linspace(0,2*pi,30); % Coordenada angular (theta)
[R,A] = meshgrid(r,a); % Mallado en (rho,theta)

% Conversión a coordenadas cartesianas
X = R.*cos(A);
Y = R.*sin(A);

% Nueva función potencial
phi = (R + 1./R).*cos(A) + A/(4*pi);

figure
surf(X,Y,phi);
shading interp
hold on
plot(1*cos(a),1*sin(a),'k','LineWidth',2); % Obstáculo
axis([-4 4 -4 4])
title('Nueva función potencial')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off

% Campo de velocidades (en cartesianas) sobre curvas de nivel de φ
figure
contour(X,Y,phi,30);
hold on

u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

quiver(X,Y,Cx,Cy)
plot(1*cos(a),1*sin(a),'k','LineWidth',1)
axis equal
axis([-4 4 -4 4])
title('Campo de velocidades (nuevo potencial)')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off
}}

Gracias a la representación anterior se observa que los vectores velocidad siguen siendo ortogonales a las líneas equipotenciales del nuevo potencial.

===Rotacional nulo y divergencia nula del nuevo campo===

Al ser <math>\vec u</math> el gradiente de un potencial escalar, su rotacional sigue siendo nulo:

<center><math>\nabla\times\vec u = \vec 0.</math></center>

En coordenadas cilíndricas se verifica explícitamente mediante la expresión

<center><math>
\nabla\times\vec u =
\frac{1}{\rho}
\begin{vmatrix}
\vec e_\rho & \vec e_\theta & \vec e_z\\
\frac{\partial}{\partial\rho} &
\frac{\partial}{\partial\theta} &
\frac{\partial}{\partial z}\\
u_\rho & \rho u_\theta & 0
\end{vmatrix}
= \vec 0,
</math></center>

lo que muestra que el flujo sigue siendo irrotacional.

Por otro lado, la divergencia del nuevo campo de velocidades viene dada por

<center><math>
\nabla\cdot\vec u(\rho,\theta,z)
=\frac{1}{\rho}\left\{
\frac{\partial}{\partial\rho}(\rho u_\rho)
+\frac{\partial}{\partial\theta}(u_\theta)
+\frac{\partial}{\partial z}(\rho u_z)
\right\},
</math></center>

y sustituyendo las expresiones de <math>u_\rho</math> y <math>u_\theta</math> se obtiene de nuevo

<center><math>\nabla\cdot\vec u = 0.</math></center>

La divergencia nula indica que el flujo sigue siendo incompresible: no hay creación ni destrucción de volumen de fluido en el dominio.

===Líneas de corriente del nuevo flujo===

Para determinar las líneas de corriente asociadas al nuevo campo de velocidades, consideramos el campo <math>\vec v</math> ortogonal a <math>\vec u</math>, definido por

<center><math>\vec v = \vec k \times \vec u.</math></center>

En coordenadas cilíndricas se tiene

<center>
<math>
\vec v
= \left(\left(1+\tfrac{1}{\rho^2}\right)\sin\theta - \tfrac{1}{4\pi\rho}\right)\vec e_\rho
+ \left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta.
</math>
</center>

Dado que el rotacional de <math>\vec v</math> también es nulo, existe un potencial escalar <math>\psi</math> tal que

<center><math>\nabla\psi = \vec v.</math></center>

A partir de las componentes de <math>\vec v</math> se obtiene el siguiente sistema:

<center><math>
\frac{\partial\psi}{\partial\rho}
= \left(1+\frac{1}{\rho^2}\right)\sin\theta - \frac{1}{4\pi\rho},\qquad
\frac{1}{\rho}\frac{\partial\psi}{\partial\theta}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta.
</math></center>

Integrando la primera ecuación respecto de <math>\rho</math> y comprobando la compatibilidad con la segunda, llegamos a

<center><math>
\psi(\rho,\theta)
= \left(\rho-\frac{1}{\rho}\right)\sin\theta
-\frac{1}{4\pi}\ln\rho + C,
</math></center>

donde <math>C</math> es una constante. Las curvas <math>\psi=\text{cte}</math> representan las líneas de corriente del nuevo flujo. Su representación numérica muestra cómo la introducción del término angular modifica el patrón de las corrientes, introduciendo una ligera asimetría respecto al caso inicial.

[[Archivo:lineasdecorrienteparaelnuevopotencial.png|thumb|400px|right|Líneas de corriente para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.3: Líneas de corriente del nuevo flujo

r = linspace(1,5,80);
a = linspace(0,2*pi,80);
[R,A] = meshgrid(r,a);

X = R.*cos(A);
Y = R.*sin(A);

% Campo de velocidades (reutilizamos las expresiones de 9.1)
u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

% Función de corriente psi
psi = (R - 1./R).*sin(A) - (1/(4*pi))*log(R);

figure
contourf(X,Y,psi,30,'LineStyle','-');
hold on
colorbar
quiver(X(1:3:end,1:3:end),Y(1:3:end,1:3:end), ...
Cx(1:3:end,1:3:end),Cy(1:3:end,1:3:end),1.3,'w');
plot(cos(a),sin(a),'k','LineWidth',2)
axis equal
axis([-4 4 -4 4])
title('Líneas de corriente para el nuevo potencial')
xlabel('EJE X'); ylabel('EJE Y');
hold off
}}

==Poster==
En el poster a continuación se muestra resumidamente el trabajo.
<center>[[Archivo:posterg68.JPG]]</center>

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-10T09:52:55Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

{{ TrabajoED | Flujo alrededor de un obstáculo circular. Grupo 68| [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] | Nicolás Rus Pérez, Jeremy García Herrera, Antonio Rodriguez Montes, Franco Luis Santillan, Juan Jose Salhua Palma.}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]
Un fluido es un medio contínuo deformable con la propiedad de que no transmite tensiones tangenciales. Estos no tienen una forma definida y por tanto adoptan la forma del recipiente en el que se encuentran, y pueden encontrarse tanto en estado líquido (agua) como gaseoso (aire).

En este trabajo se estudiará el flujo de un fluido alrededor de un obstáculo circular considerando a este un fluido incompresible. Esta propiedad indica que el fluido no disminuye su volumen al serle aplicado una fuerza, y es propia de los fluidos en estado líquido.

Por conveniencia se trabajará en coordenadas con base cilíndrica (polares).

==Mallado==
Se genera a continuación un mallado en el que se representan los puntos interiores a la región ocupada por el fluido. Por convenio se representará el obstáculo circular centrado en el origen de coordenadas y con radio 1. De esta forma la región ocupada por el fluido será todo el exterior de la circunferencia unidad.

El objetivo de realizar este mallado es dividir el medio continuo (la región ocupada por el fluido) en subdominios más pequeños denominados elementos finitos. El conjunto de estos elementos forman el mallado que cubre toda la región y permite realizar cálculos y simulaciones numéricas particulares en cada uno de ellos.

Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo <math>[-4,4]×[-4,4]</math>.
[[Archivo:Mallado1g68.JPG|thumb|right|Región del fluido y obstáculo]]
{{matlab|codigo=

r=linspace(1,5,50); %Definen el rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Genera el mallado

hold on
X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);
Z=0.*R;
mesh(X,Y,0*Z); %Dibuja la región
plot(cos(a),sin(a),'k','lineWidth',1); %Dibuja el obstáculo
axis([-4 4,-4 4]); %Se fijan los ejes
title ('Regíon Fluido'); %Definen los ejes y el nombre de la gráfica
xlabel 'EJE X'
ylabel 'EJE Y'
hold off
}}

==Campo de velocidad de las partículas==
Una vez definido el mallado de trabajo del fluido y el obstáculo, se procede a calcular las velocidades de las partículas del fluido en los distintos puntos de la región ocupada por el fluido incompresible. Partimos sabiendo que el campo de velocidades viene definido por el gradiente de la siguiente función potencial:

<center><math>\varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) </math></center>

A continuación se representa graficamente la función potencial:
[[Archivo:Mallado222g68.JPG|thumb|400px|right|Función potencial en la región del fluido]]
{{matlab|codigo=

r=linspace(1,5,30); %Define rho (r) y theta (a)
a=linspace(0,2*pi,30);

[R,A]=meshgrid(r,a); %CGenra el mallado

hold on
X=R.*cos(A); %Parametrización la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Define la función potencial
Z=f(R,A);
surf(X,Y,Z); %Representa la función
plot(1*cos(a),1*sin(a),'k','lineWidth',3); %Representa el obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de colores, define el titulo de la gráfica y los ejes
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
hold off
}}
Para calcular el campo de velocidades de las partículas del fluido es necesario el campo gradiente de la funcion potencial anterior, y dicho gradiente viene definido por la siguiente expresión:

<center>\(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)</center>

Aplicando dicha expresión a la función potencial, el gradiente resultante sera:

<center><math> \vec u(\vec e_\rho,\vec e_\theta,\vec e_z)=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta)]\vec e_\theta </math></center>

Dado que la representación en MATLAB se realiza en coordenadas cartesianas <math>(\vec i,\vec j,\vec k)</math>, y el gradiente calculado se encuentra en coordenadas ciliíndricas <math>(\vec e_\rho,\vec e_\theta,\vec e_z)</math> ,debemos de pasar el gradiente de la función potencial a las primeras. Esto es posible gracias a la matriz de cambio de base entre coordenadas cilíndricas y cartesianas, por lo que la multiplicamos por el vector gradiente dado <math>(\vec u)</math> en coordenadas cilíndricas.

Matriz de cambio de base de coordenadas cilindricas a coordenadas cartesianas:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z\\
\end{pmatrix}
</center>

Por lo que:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
(1-\frac{1}{\rho^2})\cos(\theta)\\
-(1+\frac{1}{\rho^2})\sin(\theta)\\
0\\
\end{pmatrix}
</center>
Resolviendo a coordenadas cartesianas:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta)\right) \vec j </math></center>

A continuación representamos el campo de velocidades de las partículas sobre las superficies equipotenciales del flujo del fluido.
[[Archivo:Mallado3g68.JPG|thumb|500px|right|Campo de velocidades de las partículas sobre las líneas equipotenciales del fluido]]
{{matlab|codigo=

r=linspace(1,5,50); %Define rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Crea la malla

X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Definimos la función potencial
Z=f(R,A); %Aplica la función
contour(X,Y,Z,30); %Dibuja las curvas de nivel
hold on %Define las componentes X e Y del gradiente
Cx=(1-(1./R.^2)).*cos(A).^2+(1+(1./R.^2)).*sin(A).^2;
Cy=(1-(1./R.^2)).*sin(A).*cos(A)-(1+(1./R.^2)).*sin(A).*cos(A);
quiver(X,Y,Cx,Cy); %Dibuja el campo de velocidades
plot(1*cos(a),1*sin(a),'k','lineWidth',1); %Representa del obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de color, titulo y nombre a los ejes
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}
Gracias a la gráfica anterior se observa como los vectores de velocidad definidos por <math> \vec u </math> son siempre perpendiculares a las líneas equipotenciales del fluido.
<center>[[Archivo:Mallado33.JPG]]</center>

==Rotacional Nulo y Divergencia Nula==
===Rotacional Nulo===
El rotacional nulo del campo de velocidades indica que el flujo es '''irrotacional''' , es decir , carece de vórtices o giros locales. Esto significa que una partícula de fluido no rota sobre sí misma mientras se desplaza, característica típica de los flujos potenciales.

El operador rotacional nos indica la dirección y la velocidad del giro del campo en cada punto.

A continuación se demuestra que el rotacional es nulo:

<center><math>
\nabla\times \vec{u} = \frac{1}{\rho} \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & \rho(-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix}=</math>

<math>
= \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} \left[ \rho \left( -(1 + \frac{1}{{\rho^2}})\sin(\theta) \right) \right]\vec{e}_z - \frac{\partial}{\partial \theta} \left[ \cos(\theta) \left(1 - \frac{1}{\rho^2}\right) \right] \vec{e}_z \right]=\frac{1}{\rho} \left[ -(1 - \frac{1}{\rho^2})\sin(\theta) + \sin(\theta)(1 - \frac{1}{\rho^2}) \right]=0
</math>

<math>\boxed { \nabla \times \bar { u } =0 } </math></center>

Se obtiene que el rotacional es 0, lo cual indica que las partículas del fluido no están girando

===Divergencia Nula===
Interpretación de la divergencia nula: La condición ∇·u = 0 corresponde a un flujo incompresible, donde la densidad del fluido permanece constante y el volumen de cualquier elemento material no cambia a lo largo del tiempo. Esto refleja la conservación local de la masa en ausencia de fuentes o sumideros.

<center><math>
\nabla \cdot \vec{u}(\rho, \theta, z) = \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} (\rho u_\rho) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (\rho u_z) \right\}= </math>

<math>
= \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} \left( (\rho - \frac{1}{\rho})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) \right) \right\}=
</math>
<math>
\frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0
</math>

<math>\boxed { div(\bar { u } )=0 } </math></center>

== Líneas de corriente ==
Las líneas de corriente del flujo alrededor de un cilindro circular (en sección transversal) presentan un patrón simétrico con respecto al eje horizontal (eje x) y antisimétrico con respecto al eje vertical (eje y).
[[Archivo:Linea_corrienteTangencial.png|thumb|400px|right|Líneas de corrientes]]

Las líneas de corriente <math>\psi = \text{cte}</math> describen las trayectorias del flujo potencial alrededor del obstáculo circular, mostrando simetría y adaptación suave a la geometría. Estas líneas son tangentes al campo de velocidades en cada punto, verificando que <math>\vec{u}\cdot\nabla\psi = 0</math> . El patrón obtenido revela aceleración del fluido cerca del cilindro, condición de no penetración en la superficie, y recuperación del flujo uniforme lejos del obstáculo, características típicas de un flujo incompresible e irrotacional.

<math>\vec{v}</math> se calculará con : <math>\vec{v}=\vec{k}\times \vec{u}</math>

<math>\vec{v}</math> es irrotacional dadon que el rotacional calculado en el apartado 3.1 es nulo.
De igual forma iene un potencial escalar <math>\psi </math> el cual tiene <math>\vec{v}</math> de gradiente.

Para el potencial escalar <math>\psi </math> se resolverá este sistema de ecuaciones: <math>\nabla\psi =\vec{v}</math> cumpliendo la condición anterior de líneas <math>\psi =cte</math>.

Empezamos las operaciones:
<math>\vec{v} =\vec{k}\times \vec{u} = \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
0 & 0 & 1 \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & (-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

<center><math>\boxed{\frac{\partial \psi}{\partial \rho} = V\rho,
\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = V\theta,
\frac{\partial \psi}{\partial z} = Vz}
</math></center>

<math>\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho = \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

<math>\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})\Rightarrow\psi(\rho,\theta,Z)= \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

Dando como resultado
<math>\psi: \sin{\theta}(\rho-\frac{1}{\rho})=cte</math>
[[Archivo:Psi_lineasCorriente.png|thumb|400px|right|Líneas de corrientes]]
{{matlab|codigo=
% Definir dominio en coordenadas cilíndricas
r = linspace(1, 5, 50); % Coordenada radial (ρ) desde r=1 (obstáculo)
a = linspace(0, 2*pi, 50); % Coordenada angular (θ)

% Crear malla
[R, A] = meshgrid(r, a);

% Convertir a coordenadas cartesianas
X = R .* cos(A);
Y = R .* sin(A);

% 1. FUNCIÓN POTENCIAL φ (phi)
phi = (R + (1./R)) .* cos(A);
% 2. CAMPO DE VELOCIDADES u = ∇φ

% Componentes en coordenadas cilíndricas:
% u_ρ = (1 - 1/ρ²)cos(θ)
% u_θ = -(1 + 1/ρ²)sin(θ)

u_rho = (1 - (1./R.^2)) .* cos(A);
u_theta = -(1 + (1./R.^2)) .* sin(A);

% Convertir a coordenadas cartesianas
u_x = u_rho .* cos(A) - u_theta .* sin(A);
u_y = u_rho .* sin(A) + u_theta .* cos(A);

% 3. FUNCIÓN DE CORRIENTE ψ (psi)
psi = (R - 1./R) .* sin(A);

% 4. GRÁFICAS

figure('Position', [100, 100, 1200, 800]);

% Graf 1: Líneas de corriente (ψ) y campo de velocidades
subplot(2, 2, 2);
contourf(X, Y, psi, 30, 'LineStyle', '-', 'LineWidth', 0.5);
hold on;
colorbar;
title('Líneas de Corriente (ψ) y Campo de Velocidades');
xlabel('EJE X');
ylabel('EJE Y');
axis equal;
axis([-4, 4, -4, 4]);
grid on;

% Dibujar campo de velocidades
quiver(X(1:3:end, 1:3:end), Y(1:3:end, 1:3:end), ...
u_x(1:3:end, 1:3:end), u_y(1:3:end, 1:3:end), ...
1.5, 'w', 'LineWidth', 0.7, 'MaxHeadSize', 0.5);

% Dibujar obstáculo circular
plot(cos(theta_circ), sin(theta_circ), 'k', 'LineWidth', 2, ...
'Color', [0.2, 0.2, 0.2]);
fill(cos(theta_circ), sin(theta_circ), [0.8, 0.8, 0.8], ...
'EdgeColor', [0.2, 0.2, 0.2], 'LineWidth', 1.5);
hold off;
}}

== Velocidad extrema y puntos de remanso==

En la frontera del obstáculo, definida por la circunferencia unidad <math>\rho=1</math>, el módulo de la velocidad viene dado por

<center><math>|\vec u| = 2|\sin\theta|</math></center>

De esta expresión se deduce que la velocidad máxima se alcanza cuando

<math>\sin\theta=\pm1</math>, es decir, en

<center><math>\theta=\frac{\pi}{2},\ \frac{3\pi}{2}</math></center>

que corresponden a los puntos cartesianos

<center><math>(0,1)\ \text{y}\ (0,-1)</math></center>

Estos son los puntos donde la velocidad del fluido es máxima sobre la frontera.

Por otro lado, la velocidad se anula cuando

<center><math>\sin\theta=0 \Rightarrow \theta=0,\ \pi</math></center>

lo que corresponde a los puntos

<center><math>(1,0)\ \text{y}\ (-1,0)</math></center>

Estos puntos reciben el nombre de '''puntos de remanso'''.

Para verificarlo numéricamente en MATLAB se ha representado el módulo de la velocidad sobre la frontera:
[[Archivo:Subir1.png|thumb|500px|center|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular. La velocidad es nula en los puntos de remanso y máxima en la dirección vertical.]]
{{matlab|codigo=
nr=100; nt=200;
r=linspace(1,5,nr);
th=linspace(0,2*pi,nt);
[R,TH]=meshgrid(r,th);

ur=(1-1./R.^2).*cos(TH);
uth=-(1+1./R.^2).*sin(TH);

Ux=ur.*cos(TH)-uth.*sin(TH);
Uy=ur.*sin(TH)+uth.*cos(TH);

Umag=sqrt(Ux.^2+Uy.^2);

UmagS=Umag(:,1);
THS=TH(:,1);

[Umax,imax]=max(UmagS);
[Umin,imin]=min(UmagS);

thetaMax=THS(imax);
thetaMin=THS(imin);

figure
plot(THS,UmagS), grid on
hold on
plot(thetaMax,Umax,'ro',thetaMin,Umin,'bo')
title('Modulo de la velocidad sobre la frontera S')
xlabel('\theta'), ylabel('|\itu|')
hold off
}}

== Presión del fluido y ecuación de Bernoulli==

Se supone que la densidad del fluido es constante e igual a

<center><math>\rho=2</math></center>

y se cumple la ecuación de Bernoulli

<center><math>\frac12\rho|\vec u|^2+p=C</math></center>

Despejando la presión se obtiene

<center><math>p=C-\frac12\rho|\vec u|^2</math></center>

Tomando para la constante el valor <math>C=1</math>, la presión queda completamente determinada en toda la región del fluido.

El campo de presiones se ha representado en MATLAB mediante:

{{matlab|codigo=
rho=2;
C=1;

P=C-0.5*rho*Umag.^2;

PS=P(:,1);
[Pmax,ipmax]=max(PS);
[Pmin,ipmin]=min(PS);

thetaPmax=THS(ipmax);
thetaPmin=THS(ipmin);

figure
contourf(R.*cos(TH),R.*sin(TH),P,30,'LineStyle','none');
colorbar
hold on
viscircles([0 0],1,'Color','k')
axis equal
title('Presión p(x,y)')
xlabel('x'), ylabel('y')
hold off
}}
{| class="wikitable" style="border:0; margin:auto; background:none;"
| style="border:0;" | [[Archivo:Subir3.png|thumb|500px|Campo ´de presiones según la ecuación de Bernoulli]]
| style="border:0;" | [[Archivo:Subir2.png|thumb|500px|Módulo de la velocidad ∣𝑢⃗∣ sobre la frontera del obstáculo 𝑟=1 en función del ángulo 𝜃]]
|}
Se observa que la presión es máxima en los puntos de remanso, donde la velocidad es nula, y mínima en los puntos donde la velocidad alcanza su valor máximo. Esto concuerda con la ecuación de Bernoulli, ya que un aumento del módulo de la velocidad implica necesariamente una disminución de la presión.

==Velocidad y presión al rodear el obstáculo.==

Ya teniendo el campo de presiones se puede apreciar que tiene relación con la velocidad, ya que cuando disminuye la presión aumenta el módulo de la velocidad, tratándose de un resultado lógico debido a que cuando hablamos de un fluido incompresible su valor de Trinomio de Bernoulli ha de ser constante.
Lo podemos ver mediante la siguiente expresión: (Ec. Bernoulli).

<math>B=z+P/γ+v^2/2g</math>

B=valor del trinomio de Bernoulli

P= presión

z= altura

v=velocidad que lleva el fluido en determinado punto

g=valor de la gravedad

γ=valor del peso específico

Donde podemos ver que si la presión aumenta la velocidad tiene que disminuir y viceversa.

==Paradoja de D'Alembert.==

La magnitud pn representa la fuerza que el fluido ejerce de forma puntual sobre la superficie del obstáculo. Para estudiar si existe empuje neto en la dirección horizontal, basta analizar la suma de las proyecciones de todas esas fuerzas sobre el vector i. Si dicha suma resulta nula, concluiremos que no hay fuerza total en esa dirección.

Del apartado anterior conocemos ya la presión en la frontera del círculo unidad:
<math>p=10-∥u∥^2=9-4 sin⁡^2 θ-4sin⁡θ</math>

Cálculo del vector normal

Sobre la circunferencia de radio 1 el vector normal exterior es

<math>n=-eρ.</math>

Para expresarlo en coordenadas cartesianas utilizamos la matriz de cambio de base de polares a cartesianas. Tomando su traspuesta se obtiene:

\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}

Aplicando esta matriz a <math>(1,0,0)^T</math>, que corresponde al versor radial eρ, se obtiene:
\begin{pmatrix}
\cos \theta &\\
sen \theta &\\
0
\end{pmatrix}
<math>=cos⁡θi+sin⁡θj</math>

Al tomar el signo negativo para el normal exterior,

<math>n=cos⁡θi-sin⁡θj</math>,

de donde su proyección sobre i es:

<math>n⋅i=-cos⁡θ</math>.

Integración sobre la frontera

El borde del dominio se parametriza como
<math>(ρ,θ)=(1,θ), con θ∈[0,2π)</math>.

La contribución total de la presión en la dirección horizontal es:
<math>\int_{0}^{2π}p(n⋅i)dθ=\int_{0}^{2π}p(-cos⁡θ)dθ.</math>

Sustituyendo la expresión de p:
<math>\int_{0}^{2π}(9-4sin⁡^2 θ-4sin⁡θ)(-cos⁡θ)dθ.</math>

Integramos término a término:
<math>\int_{0}^{2π}-9cos⁡θdθ+\int_{0}^{2π}4sin⁡^2 θcos⁡θdθ+\int_{0}^{2π}4sin⁡θcos⁡θdθ.</math>

Cada una de estas integrales es periódica con una primitiva continua en todo el intervalo; evaluándolas entre 0 y 2π se obtiene:
<math>[-9sin⁡θ]_{0}^{2π}=0,[(4sin⁡^3 θ)/3]_{0}^{2π}=0,[2sin⁡^2 θ]_{0}^{2π}=0.</math>

Por tanto,
<math>\int_{0}^{2π}p(n⋅i)dθ=0.</math>

Conclusión

La componente horizontal de la fuerza total ejercida por el fluido sobre el obstáculo es nula. En consecuencia, el fluido no produce empuje neto en la dirección i.

==Cambio de función potencial==

En este apartado repetimos los procesos de los puntos (2), (3) y (4), pero tomando ahora como función potencial

<center><math>\varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}</math></center>

El obstáculo sigue siendo la circunferencia unidad centrada en el origen, y el fluido ocupa la región exterior, trabajándose de nuevo en coordenadas cilíndricas (polares).

===Representación de la función potencial y del campo de velocidades===

En primer lugar representamos la nueva función potencial <math>\varphi</math> sobre la región ocupada por el fluido, manteniendo el mismo mallado del apartado inicial. Esta representación permite visualizar cómo el término angular <math>\theta/(4\pi)</math> introduce una ligera inclinación en la superficie potencial respecto al caso anterior.

[[Archivo:nuevafuncionpotencial.png|thumb|400px|right|Nueva función potencial del flujo con término angular]]

La velocidad de las partículas del fluido viene dada de nuevo por el gradiente del potencial. En coordenadas cilíndricas, el campo de velocidades es

<center><math>\vec u(\rho,\theta,z)=\nabla\varphi
= u_\rho\,\vec e_\rho+u_\theta\,\vec e_\theta+u_z\,\vec e_z</math></center>

donde

<center><math>
u_\rho=\frac{\partial\varphi}{\partial\rho}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta, \qquad
u_\theta=\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}
= -\left(1+\frac{1}{\rho^2}\right)\sin\theta+\frac{1}{4\pi\rho},\qquad
u_z=0.
</math></center>

Para poder representarlo en MATLAB se pasa el campo a coordenadas cartesianas mediante la matriz de cambio de base de cilíndricas a cartesianas:

<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k
\end{pmatrix}
=
\begin{pmatrix}
\cos\theta & -\sin\theta & 0\\
\sin\theta & \cos\theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z
\end{pmatrix}
</center>

Aplicando esta matriz al vector <math>(u_\rho,u_\theta,0)^T</math> se obtiene

<center><math>
\vec u(\vec i,\vec j,\vec k)
=\left[(1-\tfrac{1}{\rho^2})\cos^2\theta+(1+\tfrac{1}{\rho^2})\sin^2\theta
-\frac{\sin\theta}{4\pi\rho}\right]\vec i
+\left[(1-\tfrac{1}{\rho^2})\sin\theta\cos\theta
-(1+\tfrac{1}{\rho^2})\sin\theta\cos\theta
+\frac{\cos\theta}{4\pi\rho}\right]\vec j.
</math></center>

Sobre las curvas de nivel de <math>\varphi</math> se representa este nuevo campo de velocidades:

[[Archivo:campodevelocidadesnuevopotencial.png|thumb|500px|right|Campo de velocidades sobre las líneas equipotenciales para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.1: Potencial y campo de velocidades para la nueva función φ

r = linspace(1,5,30); % Coordenada radial (rho)
a = linspace(0,2*pi,30); % Coordenada angular (theta)
[R,A] = meshgrid(r,a); % Mallado en (rho,theta)

% Conversión a coordenadas cartesianas
X = R.*cos(A);
Y = R.*sin(A);

% Nueva función potencial
phi = (R + 1./R).*cos(A) + A/(4*pi);

figure
surf(X,Y,phi);
shading interp
hold on
plot(1*cos(a),1*sin(a),'k','LineWidth',2); % Obstáculo
axis([-4 4 -4 4])
title('Nueva función potencial')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off

% Campo de velocidades (en cartesianas) sobre curvas de nivel de φ
figure
contour(X,Y,phi,30);
hold on

u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

quiver(X,Y,Cx,Cy)
plot(1*cos(a),1*sin(a),'k','LineWidth',1)
axis equal
axis([-4 4 -4 4])
title('Campo de velocidades (nuevo potencial)')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off
}}

Gracias a la representación anterior se observa que los vectores velocidad siguen siendo ortogonales a las líneas equipotenciales del nuevo potencial.

===Rotacional nulo y divergencia nula del nuevo campo===

Al ser <math>\vec u</math> el gradiente de un potencial escalar, su rotacional sigue siendo nulo:

<center><math>\nabla\times\vec u = \vec 0.</math></center>

En coordenadas cilíndricas se verifica explícitamente mediante la expresión

<center><math>
\nabla\times\vec u =
\frac{1}{\rho}
\begin{vmatrix}
\vec e_\rho & \vec e_\theta & \vec e_z\\
\frac{\partial}{\partial\rho} &
\frac{\partial}{\partial\theta} &
\frac{\partial}{\partial z}\\
u_\rho & \rho u_\theta & 0
\end{vmatrix}
= \vec 0,
</math></center>

lo que muestra que el flujo sigue siendo irrotacional.

Por otro lado, la divergencia del nuevo campo de velocidades viene dada por

<center><math>
\nabla\cdot\vec u(\rho,\theta,z)
=\frac{1}{\rho}\left\{
\frac{\partial}{\partial\rho}(\rho u_\rho)
+\frac{\partial}{\partial\theta}(u_\theta)
+\frac{\partial}{\partial z}(\rho u_z)
\right\},
</math></center>

y sustituyendo las expresiones de <math>u_\rho</math> y <math>u_\theta</math> se obtiene de nuevo

<center><math>\nabla\cdot\vec u = 0.</math></center>

La divergencia nula indica que el flujo sigue siendo incompresible: no hay creación ni destrucción de volumen de fluido en el dominio.

===Líneas de corriente del nuevo flujo===

Para determinar las líneas de corriente asociadas al nuevo campo de velocidades, consideramos el campo <math>\vec v</math> ortogonal a <math>\vec u</math>, definido por

<center><math>\vec v = \vec k \times \vec u.</math></center>

En coordenadas cilíndricas se tiene

<center>
<math>
\vec v
= \left(\left(1+\tfrac{1}{\rho^2}\right)\sin\theta - \tfrac{1}{4\pi\rho}\right)\vec e_\rho
+ \left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta.
</math>
</center>

Dado que el rotacional de <math>\vec v</math> también es nulo, existe un potencial escalar <math>\psi</math> tal que

<center><math>\nabla\psi = \vec v.</math></center>

A partir de las componentes de <math>\vec v</math> se obtiene el siguiente sistema:

<center><math>
\frac{\partial\psi}{\partial\rho}
= \left(1+\frac{1}{\rho^2}\right)\sin\theta - \frac{1}{4\pi\rho},\qquad
\frac{1}{\rho}\frac{\partial\psi}{\partial\theta}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta.
</math></center>

Integrando la primera ecuación respecto de <math>\rho</math> y comprobando la compatibilidad con la segunda, llegamos a

<center><math>
\psi(\rho,\theta)
= \left(\rho-\frac{1}{\rho}\right)\sin\theta
-\frac{1}{4\pi}\ln\rho + C,
</math></center>

donde <math>C</math> es una constante. Las curvas <math>\psi=\text{cte}</math> representan las líneas de corriente del nuevo flujo. Su representación numérica muestra cómo la introducción del término angular modifica el patrón de las corrientes, introduciendo una ligera asimetría respecto al caso inicial.

[[Archivo:lineasdecorrienteparaelnuevopotencial.png|thumb|400px|right|Líneas de corriente para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.3: Líneas de corriente del nuevo flujo

r = linspace(1,5,80);
a = linspace(0,2*pi,80);
[R,A] = meshgrid(r,a);

X = R.*cos(A);
Y = R.*sin(A);

% Campo de velocidades (reutilizamos las expresiones de 9.1)
u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

% Función de corriente psi
psi = (R - 1./R).*sin(A) - (1/(4*pi))*log(R);

figure
contourf(X,Y,psi,30,'LineStyle','-');
hold on
colorbar
quiver(X(1:3:end,1:3:end),Y(1:3:end,1:3:end), ...
Cx(1:3:end,1:3:end),Cy(1:3:end,1:3:end),1.3,'w');
plot(cos(a),sin(a),'k','LineWidth',2)
axis equal
axis([-4 4 -4 4])
title('Líneas de corriente para el nuevo potencial')
xlabel('EJE X'); ylabel('EJE Y');
hold off
}}

==Poster==
En el poster a continuación se muestra resumidamente el trabajo.
<center>[[Archivo:posterg68.JPG]]</center>

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-07T23:03:26Z

Jeremy.garcia.herrera: /* Líneas de corriente del nuevo flujo */

{{ TrabajoED | Flujo alrededor de un obstáculo circular. Grupo 68| [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] | Nicolás Rus Pérez, Jeremy García Herrera, Antonio Rodriguez Montes, Franco Luis Santillan, Juan Jose Salhua Palma.}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]
Un fluido es un medio contínuo deformable con la propiedad de que no transmite tensiones tangenciales. Estos no tienen una forma definida y por tanto adoptan la forma del recipiente en el que se encuentran, y pueden encontrarse tanto en estado líquido (agua) como gaseoso (aire).

En este trabajo se estudiará el flujo de un fluido alrededor de un obstáculo circular considerando a este un fluido incompresible. Esta propiedad indica que el fluido no disminuye su volumen al serle aplicado una fuerza, y es propia de los fluidos en estado líquido.

Por conveniencia se trabajará en coordenadas con base cilíndrica (polares).

==Mallado==
Se genera a continuación un mallado en el que se representan los puntos interiores a la región ocupada por el fluido. Por convenio se representará el obstáculo circular centrado en el origen de coordenadas y con radio 1. De esta forma la región ocupada por el fluido será todo el exterior de la circunferencia unidad.

El objetivo de realizar este mallado es dividir el medio continuo (la región ocupada por el fluido) en subdominios más pequeños denominados elementos finitos. El conjunto de estos elementos forman el mallado que cubre toda la región y permite realizar cálculos y simulaciones numéricas particulares en cada uno de ellos.

Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo <math>[-4,4]×[-4,4]</math>.
[[Archivo:Mallado1g68.JPG|thumb|right|Región del fluido y obstáculo]]
{{matlab|codigo=

r=linspace(1,5,50); %Definen el rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Genera el mallado

hold on
X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);
Z=0.*R;
mesh(X,Y,0*Z); %Dibuja la región
plot(cos(a),sin(a),'k','lineWidth',1); %Dibuja el obstáculo
axis([-4 4,-4 4]); %Se fijan los ejes
title ('Regíon Fluido'); %Definen los ejes y el nombre de la gráfica
xlabel 'EJE X'
ylabel 'EJE Y'
hold off
}}

==Campo de velocidad de las partículas==
Una vez definido el mallado de trabajo del fluido y el obstáculo, se procede a calcular las velocidades de las partículas del fluido en los distintos puntos de la región ocupada por el fluido incompresible. Partimos sabiendo que el campo de velocidades viene definido por el gradiente de la siguiente función potencial:

<center><math>\varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) </math></center>

A continuación se representa graficamente la función potencial:
[[Archivo:Mallado222g68.JPG|thumb|400px|right|Función potencial en la región del fluido]]
{{matlab|codigo=

r=linspace(1,5,30); %Define rho (r) y theta (a)
a=linspace(0,2*pi,30);

[R,A]=meshgrid(r,a); %CGenra el mallado

hold on
X=R.*cos(A); %Parametrización la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Define la función potencial
Z=f(R,A);
surf(X,Y,Z); %Representa la función
plot(1*cos(a),1*sin(a),'k','lineWidth',3); %Representa el obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de colores, define el titulo de la gráfica y los ejes
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
hold off
}}
Para calcular el campo de velocidades de las partículas del fluido es necesario el campo gradiente de la funcion potencial anterior, y dicho gradiente viene definido por la siguiente expresión:

<center>\(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)</center>

Aplicando dicha expresión a la función potencial, el gradiente resultante sera:

<center><math> \vec u(\vec e_\rho,\vec e_\theta,\vec e_z)=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta)]\vec e_\theta </math></center>

Dado que la representación en MATLAB se realiza en coordenadas cartesianas <math>(\vec i,\vec j,\vec k)</math>, y el gradiente calculado se encuentra en coordenadas ciliíndricas <math>(\vec e_\rho,\vec e_\theta,\vec e_z)</math> ,debemos de pasar el gradiente de la función potencial a las primeras. Esto es posible gracias a la matriz de cambio de base entre coordenadas cilíndricas y cartesianas, por lo que la multiplicamos por el vector gradiente dado <math>(\vec u)</math> en coordenadas cilíndricas.

Matriz de cambio de base de coordenadas cilindricas a coordenadas cartesianas:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z\\
\end{pmatrix}
</center>

Por lo que:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
(1-\frac{1}{\rho^2})\cos(\theta)\\
-(1+\frac{1}{\rho^2})\sin(\theta)\\
0\\
\end{pmatrix}
</center>
Resolviendo a coordenadas cartesianas:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta)\right) \vec j </math></center>

A continuación representamos el campo de velocidades de las partículas sobre las superficies equipotenciales del flujo del fluido.
[[Archivo:Mallado3g68.JPG|thumb|500px|right|Campo de velocidades de las partículas sobre las líneas equipotenciales del fluido]]
{{matlab|codigo=

r=linspace(1,5,50); %Define rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Crea la malla

X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Definimos la función potencial
Z=f(R,A); %Aplica la función
contour(X,Y,Z,30); %Dibuja las curvas de nivel
hold on %Define las componentes X e Y del gradiente
Cx=(1-(1./R.^2)).*cos(A).^2+(1+(1./R.^2)).*sin(A).^2;
Cy=(1-(1./R.^2)).*sin(A).*cos(A)-(1+(1./R.^2)).*sin(A).*cos(A);
quiver(X,Y,Cx,Cy); %Dibuja el campo de velocidades
plot(1*cos(a),1*sin(a),'k','lineWidth',1); %Representa del obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de color, titulo y nombre a los ejes
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}
Gracias a la gráfica anterior se observa como los vectores de velocidad definidos por <math> \vec u </math> son siempre perpendiculares a las líneas equipotenciales del fluido.
<center>[[Archivo:Mallado33.JPG]]</center>

==Rotacional Nulo y Divergencia Nula==
===Rotacional Nulo===
El rotacional nulo del campo de velocidades indica que el flujo es '''irrotacional''' , es decir , carece de vórtices o giros locales. Esto significa que una partícula de fluido no rota sobre sí misma mientras se desplaza, característica típica de los flujos potenciales.

El operador rotacional nos indica la dirección y la velocidad del giro del campo en cada punto.

A continuación se demuestra que el rotacional es nulo:

<center><math>
\nabla\times \vec{u} = \frac{1}{\rho} \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & \rho(-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix}=</math>

<math>
= \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} \left[ \rho \left( -(1 + \frac{1}{{\rho^2}})\sin(\theta) \right) \right]\vec{e}_z - \frac{\partial}{\partial \theta} \left[ \cos(\theta) \left(1 - \frac{1}{\rho^2}\right) \right] \vec{e}_z \right]=\frac{1}{\rho} \left[ -(1 - \frac{1}{\rho^2})\sin(\theta) + \sin(\theta)(1 - \frac{1}{\rho^2}) \right]=0
</math>

<math>\boxed { \nabla \times \bar { u } =0 } </math></center>

Se obtiene que el rotacional es 0, lo cual indica que las partículas del fluido no están girando

===Divergencia Nula===
Interpretación de la divergencia nula: La condición ∇·u = 0 corresponde a un flujo incompresible, donde la densidad del fluido permanece constante y el volumen de cualquier elemento material no cambia a lo largo del tiempo. Esto refleja la conservación local de la masa en ausencia de fuentes o sumideros.

<center><math>
\nabla \cdot \vec{u}(\rho, \theta, z) = \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} (\rho u_\rho) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (\rho u_z) \right\}= </math>

<math>
= \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} \left( (\rho - \frac{1}{\rho})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) \right) \right\}=
</math>
<math>
\frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0
</math>

<math>\boxed { div(\bar { u } )=0 } </math></center>

== Líneas de corriente ==
Las líneas de corriente del flujo alrededor de un cilindro circular (en sección transversal) presentan un patrón simétrico con respecto al eje horizontal (eje x) y antisimétrico con respecto al eje vertical (eje y).
[[Archivo:Linea_corrienteTangencial.png|thumb|400px|right|Líneas de corrientes]]

Las líneas de corriente <math>\psi = \text{cte}</math> describen las trayectorias del flujo potencial alrededor del obstáculo circular, mostrando simetría y adaptación suave a la geometría. Estas líneas son tangentes al campo de velocidades en cada punto, verificando que <math>\vec{u}\cdot\nabla\psi = 0</math> . El patrón obtenido revela aceleración del fluido cerca del cilindro, condición de no penetración en la superficie, y recuperación del flujo uniforme lejos del obstáculo, características típicas de un flujo incompresible e irrotacional.

<math>\vec{v}</math> se calculará con : <math>\vec{v}=\vec{k}\times \vec{u}</math>

<math>\vec{v}</math> es irrotacional dadon que el rotacional calculado en el apartado 3.1 es nulo.
De igual forma iene un potencial escalar <math>\psi </math> el cual tiene <math>\vec{v}</math> de gradiente.

Para el potencial escalar <math>\psi </math> se resolverá este sistema de ecuaciones: <math>\nabla\psi =\vec{v}</math> cumpliendo la condición anterior de líneas <math>\psi =cte</math>.

Empezamos las operaciones:
<math>\vec{v} =\vec{k}\times \vec{u} = \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
0 & 0 & 1 \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & (-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

<center><math>\boxed{\frac{\partial \psi}{\partial \rho} = V\rho,
\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = V\theta,
\frac{\partial \psi}{\partial z} = Vz}
</math></center>

<math>\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho = \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

<math>\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})\Rightarrow\psi(\rho,\theta,Z)= \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

Dando como resultado
<math>\psi: \sin{\theta}(\rho-\frac{1}{\rho})=cte</math>
[[Archivo:Psi_lineasCorriente.png|thumb|400px|right|Líneas de corrientes]]
{{matlab|codigo=
% Definir dominio en coordenadas cilíndricas
r = linspace(1, 5, 50); % Coordenada radial (ρ) desde r=1 (obstáculo)
a = linspace(0, 2*pi, 50); % Coordenada angular (θ)

% Crear malla
[R, A] = meshgrid(r, a);

% Convertir a coordenadas cartesianas
X = R .* cos(A);
Y = R .* sin(A);

% 1. FUNCIÓN POTENCIAL φ (phi)
phi = (R + (1./R)) .* cos(A);
% 2. CAMPO DE VELOCIDADES u = ∇φ

% Componentes en coordenadas cilíndricas:
% u_ρ = (1 - 1/ρ²)cos(θ)
% u_θ = -(1 + 1/ρ²)sin(θ)

u_rho = (1 - (1./R.^2)) .* cos(A);
u_theta = -(1 + (1./R.^2)) .* sin(A);

% Convertir a coordenadas cartesianas
u_x = u_rho .* cos(A) - u_theta .* sin(A);
u_y = u_rho .* sin(A) + u_theta .* cos(A);

% 3. FUNCIÓN DE CORRIENTE ψ (psi)
psi = (R - 1./R) .* sin(A);

% 4. GRÁFICAS

figure('Position', [100, 100, 1200, 800]);

% Graf 1: Líneas de corriente (ψ) y campo de velocidades
subplot(2, 2, 2);
contourf(X, Y, psi, 30, 'LineStyle', '-', 'LineWidth', 0.5);
hold on;
colorbar;
title('Líneas de Corriente (ψ) y Campo de Velocidades');
xlabel('EJE X');
ylabel('EJE Y');
axis equal;
axis([-4, 4, -4, 4]);
grid on;

% Dibujar campo de velocidades
quiver(X(1:3:end, 1:3:end), Y(1:3:end, 1:3:end), ...
u_x(1:3:end, 1:3:end), u_y(1:3:end, 1:3:end), ...
1.5, 'w', 'LineWidth', 0.7, 'MaxHeadSize', 0.5);

% Dibujar obstáculo circular
plot(cos(theta_circ), sin(theta_circ), 'k', 'LineWidth', 2, ...
'Color', [0.2, 0.2, 0.2]);
fill(cos(theta_circ), sin(theta_circ), [0.8, 0.8, 0.8], ...
'EdgeColor', [0.2, 0.2, 0.2], 'LineWidth', 1.5);
hold off;
}}

== Velocidad extrema y puntos de remanso==

En la frontera del obstáculo, definida por la circunferencia unidad <math>\rho=1</math>, el módulo de la velocidad viene dado por

<center><math>|\vec u| = 2|\sin\theta|</math></center>

De esta expresión se deduce que la velocidad máxima se alcanza cuando

<math>\sin\theta=\pm1</math>, es decir, en

<center><math>\theta=\frac{\pi}{2},\ \frac{3\pi}{2}</math></center>

que corresponden a los puntos cartesianos

<center><math>(0,1)\ \text{y}\ (0,-1)</math></center>

Estos son los puntos donde la velocidad del fluido es máxima sobre la frontera.

Por otro lado, la velocidad se anula cuando

<center><math>\sin\theta=0 \Rightarrow \theta=0,\ \pi</math></center>

lo que corresponde a los puntos

<center><math>(1,0)\ \text{y}\ (-1,0)</math></center>

Estos puntos reciben el nombre de '''puntos de remanso'''.

Para verificarlo numéricamente en MATLAB se ha representado el módulo de la velocidad sobre la frontera:
[[Archivo:Subir1.png|thumb|500px|center|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular. La velocidad es nula en los puntos de remanso y máxima en la dirección vertical.]]
{{matlab|codigo=
nr=100; nt=200;
r=linspace(1,5,nr);
th=linspace(0,2*pi,nt);
[R,TH]=meshgrid(r,th);

ur=(1-1./R.^2).*cos(TH);
uth=-(1+1./R.^2).*sin(TH);

Ux=ur.*cos(TH)-uth.*sin(TH);
Uy=ur.*sin(TH)+uth.*cos(TH);

Umag=sqrt(Ux.^2+Uy.^2);

UmagS=Umag(:,1);
THS=TH(:,1);

[Umax,imax]=max(UmagS);
[Umin,imin]=min(UmagS);

thetaMax=THS(imax);
thetaMin=THS(imin);

figure
plot(THS,UmagS), grid on
hold on
plot(thetaMax,Umax,'ro',thetaMin,Umin,'bo')
title('Modulo de la velocidad sobre la frontera S')
xlabel('\theta'), ylabel('|\itu|')
hold off
}}

== Presión del fluido y ecuación de Bernoulli==

Se supone que la densidad del fluido es constante e igual a

<center><math>\rho=2</math></center>

y se cumple la ecuación de Bernoulli

<center><math>\frac12\rho|\vec u|^2+p=C</math></center>

Despejando la presión se obtiene

<center><math>p=C-\frac12\rho|\vec u|^2</math></center>

Tomando para la constante el valor <math>C=1</math>, la presión queda completamente determinada en toda la región del fluido.

El campo de presiones se ha representado en MATLAB mediante:

{{matlab|codigo=
rho=2;
C=1;

P=C-0.5*rho*Umag.^2;

PS=P(:,1);
[Pmax,ipmax]=max(PS);
[Pmin,ipmin]=min(PS);

thetaPmax=THS(ipmax);
thetaPmin=THS(ipmin);

figure
contourf(R.*cos(TH),R.*sin(TH),P,30,'LineStyle','none');
colorbar
hold on
viscircles([0 0],1,'Color','k')
axis equal
title('Presión p(x,y)')
xlabel('x'), ylabel('y')
hold off
}}
{| class="wikitable" style="border:0; margin:auto; background:none;"
| style="border:0;" | [[Archivo:Subir3.png|thumb|500px|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular]]
| style="border:0;" | [[Archivo:Subir2.png|thumb|500px|Módulo de la velocidad ∣𝑢⃗∣ sobre la frontera del obstáculo 𝑟=1 en función del ángulo 𝜃]]
|}
Se observa que la presión es máxima en los puntos de remanso, donde la velocidad es nula, y mínima en los puntos donde la velocidad alcanza su valor máximo. Esto concuerda con la ecuación de Bernoulli, ya que un aumento del módulo de la velocidad implica necesariamente una disminución de la presión.

==Velocidad y presión al rodear el obstáculo.==

Ya teniendo el campo de presiones se puede apreciar que tiene relación con la velocidad, ya que cuando disminuye la presión aumenta el módulo de la velocidad, tratándose de un resultado lógico debido a que cuando hablamos de un fluido incompresible su valor de Trinomio de Bernoulli ha de ser constante.
Lo podemos ver mediante la siguiente expresión: (Ec. Bernoulli).

<math>B=z+P/γ+v^2/2g</math>

B=valor del trinomio de Bernoulli

P= presión

z= altura

v=velocidad que lleva el fluido en determinado punto

g=valor de la gravedad

γ=valor del peso específico

Donde podemos ver que si la presión aumenta la velocidad tiene que disminuir y viceversa.

==Paradoja de D'Alembert.==

La magnitud pn representa la fuerza que el fluido ejerce de forma puntual sobre la superficie del obstáculo. Para estudiar si existe empuje neto en la dirección horizontal, basta analizar la suma de las proyecciones de todas esas fuerzas sobre el vector i. Si dicha suma resulta nula, concluiremos que no hay fuerza total en esa dirección.

Del apartado anterior conocemos ya la presión en la frontera del círculo unidad:
<math>p=10-∥u∥^2=9-4 sin⁡^2 θ-4sin⁡θ</math>

Cálculo del vector normal

Sobre la circunferencia de radio 1 el vector normal exterior es

<math>n=-eρ.</math>

Para expresarlo en coordenadas cartesianas utilizamos la matriz de cambio de base de polares a cartesianas. Tomando su traspuesta se obtiene:

\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}

Aplicando esta matriz a <math>(1,0,0)^T</math>, que corresponde al versor radial eρ, se obtiene:
\begin{pmatrix}
\cos \theta &\\
sen \theta &\\
0
\end{pmatrix}
<math>=cos⁡θi+sin⁡θj</math>

Al tomar el signo negativo para el normal exterior,

<math>n=cos⁡θi-sin⁡θj</math>,

de donde su proyección sobre i es:

<math>n⋅i=-cos⁡θ</math>.

Integración sobre la frontera

El borde del dominio se parametriza como
<math>(ρ,θ)=(1,θ), con θ∈[0,2π)</math>.

La contribución total de la presión en la dirección horizontal es:
<math>\int_{0}^{2π}p(n⋅i)dθ=\int_{0}^{2π}p(-cos⁡θ)dθ.</math>

Sustituyendo la expresión de p:
<math>\int_{0}^{2π}(9-4sin⁡^2 θ-4sin⁡θ)(-cos⁡θ)dθ.</math>

Integramos término a término:
<math>\int_{0}^{2π}-9cos⁡θdθ+\int_{0}^{2π}4sin⁡^2 θcos⁡θdθ+\int_{0}^{2π}4sin⁡θcos⁡θdθ.</math>

Cada una de estas integrales es periódica con una primitiva continua en todo el intervalo; evaluándolas entre 0 y 2π se obtiene:
<math>[-9sin⁡θ]_{0}^{2π}=0,[(4sin⁡^3 θ)/3]_{0}^{2π}=0,[2sin⁡^2 θ]_{0}^{2π}=0.</math>

Por tanto,
<math>\int_{0}^{2π}p(n⋅i)dθ=0.</math>

Conclusión

La componente horizontal de la fuerza total ejercida por el fluido sobre el obstáculo es nula. En consecuencia, el fluido no produce empuje neto en la dirección i.

==Cambio de función potencial==

En este apartado repetimos los procesos de los puntos (2), (3) y (4), pero tomando ahora como función potencial

<center><math>\varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}</math></center>

El obstáculo sigue siendo la circunferencia unidad centrada en el origen, y el fluido ocupa la región exterior, trabajándose de nuevo en coordenadas cilíndricas (polares).

===Representación de la función potencial y del campo de velocidades===

En primer lugar representamos la nueva función potencial <math>\varphi</math> sobre la región ocupada por el fluido, manteniendo el mismo mallado del apartado inicial. Esta representación permite visualizar cómo el término angular <math>\theta/(4\pi)</math> introduce una ligera inclinación en la superficie potencial respecto al caso anterior.

[[Archivo:nuevafuncionpotencial.png|thumb|400px|right|Nueva función potencial del flujo con término angular]]

La velocidad de las partículas del fluido viene dada de nuevo por el gradiente del potencial. En coordenadas cilíndricas, el campo de velocidades es

<center><math>\vec u(\rho,\theta,z)=\nabla\varphi
= u_\rho\,\vec e_\rho+u_\theta\,\vec e_\theta+u_z\,\vec e_z</math></center>

donde

<center><math>
u_\rho=\frac{\partial\varphi}{\partial\rho}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta, \qquad
u_\theta=\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}
= -\left(1+\frac{1}{\rho^2}\right)\sin\theta+\frac{1}{4\pi\rho},\qquad
u_z=0.
</math></center>

Para poder representarlo en MATLAB se pasa el campo a coordenadas cartesianas mediante la matriz de cambio de base de cilíndricas a cartesianas:

<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k
\end{pmatrix}
=
\begin{pmatrix}
\cos\theta & -\sin\theta & 0\\
\sin\theta & \cos\theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z
\end{pmatrix}
</center>

Aplicando esta matriz al vector <math>(u_\rho,u_\theta,0)^T</math> se obtiene

<center><math>
\vec u(\vec i,\vec j,\vec k)
=\left[(1-\tfrac{1}{\rho^2})\cos^2\theta+(1+\tfrac{1}{\rho^2})\sin^2\theta
-\frac{\sin\theta}{4\pi\rho}\right]\vec i
+\left[(1-\tfrac{1}{\rho^2})\sin\theta\cos\theta
-(1+\tfrac{1}{\rho^2})\sin\theta\cos\theta
+\frac{\cos\theta}{4\pi\rho}\right]\vec j.
</math></center>

Sobre las curvas de nivel de <math>\varphi</math> se representa este nuevo campo de velocidades:

[[Archivo:campodevelocidadesnuevopotencial.png|thumb|500px|right|Campo de velocidades sobre las líneas equipotenciales para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.1: Potencial y campo de velocidades para la nueva función φ

r = linspace(1,5,30); % Coordenada radial (rho)
a = linspace(0,2*pi,30); % Coordenada angular (theta)
[R,A] = meshgrid(r,a); % Mallado en (rho,theta)

% Conversión a coordenadas cartesianas
X = R.*cos(A);
Y = R.*sin(A);

% Nueva función potencial
phi = (R + 1./R).*cos(A) + A/(4*pi);

figure
surf(X,Y,phi);
shading interp
hold on
plot(1*cos(a),1*sin(a),'k','LineWidth',2); % Obstáculo
axis([-4 4 -4 4])
title('Nueva función potencial')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off

% Campo de velocidades (en cartesianas) sobre curvas de nivel de φ
figure
contour(X,Y,phi,30);
hold on

u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

quiver(X,Y,Cx,Cy)
plot(1*cos(a),1*sin(a),'k','LineWidth',1)
axis equal
axis([-4 4 -4 4])
title('Campo de velocidades (nuevo potencial)')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off
}}

Gracias a la representación anterior se observa que los vectores velocidad siguen siendo ortogonales a las líneas equipotenciales del nuevo potencial.

===Rotacional nulo y divergencia nula del nuevo campo===

Al ser <math>\vec u</math> el gradiente de un potencial escalar, su rotacional sigue siendo nulo:

<center><math>\nabla\times\vec u = \vec 0.</math></center>

En coordenadas cilíndricas se verifica explícitamente mediante la expresión

<center><math>
\nabla\times\vec u =
\frac{1}{\rho}
\begin{vmatrix}
\vec e_\rho & \vec e_\theta & \vec e_z\\
\frac{\partial}{\partial\rho} &
\frac{\partial}{\partial\theta} &
\frac{\partial}{\partial z}\\
u_\rho & \rho u_\theta & 0
\end{vmatrix}
= \vec 0,
</math></center>

lo que muestra que el flujo sigue siendo irrotacional.

Por otro lado, la divergencia del nuevo campo de velocidades viene dada por

<center><math>
\nabla\cdot\vec u(\rho,\theta,z)
=\frac{1}{\rho}\left\{
\frac{\partial}{\partial\rho}(\rho u_\rho)
+\frac{\partial}{\partial\theta}(u_\theta)
+\frac{\partial}{\partial z}(\rho u_z)
\right\},
</math></center>

y sustituyendo las expresiones de <math>u_\rho</math> y <math>u_\theta</math> se obtiene de nuevo

<center><math>\nabla\cdot\vec u = 0.</math></center>

La divergencia nula indica que el flujo sigue siendo incompresible: no hay creación ni destrucción de volumen de fluido en el dominio.

===Líneas de corriente del nuevo flujo===

Para determinar las líneas de corriente asociadas al nuevo campo de velocidades, consideramos el campo <math>\vec v</math> ortogonal a <math>\vec u</math>, definido por

<center><math>\vec v = \vec k \times \vec u.</math></center>

En coordenadas cilíndricas se tiene

<center>
<math>
\vec v
= \left(\left(1+\tfrac{1}{\rho^2}\right)\sin\theta - \tfrac{1}{4\pi\rho}\right)\vec e_\rho
+ \left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta.
</math>
</center>

Dado que el rotacional de <math>\vec v</math> también es nulo, existe un potencial escalar <math>\psi</math> tal que

<center><math>\nabla\psi = \vec v.</math></center>

A partir de las componentes de <math>\vec v</math> se obtiene el siguiente sistema:

<center><math>
\frac{\partial\psi}{\partial\rho}
= \left(1+\frac{1}{\rho^2}\right)\sin\theta - \frac{1}{4\pi\rho},\qquad
\frac{1}{\rho}\frac{\partial\psi}{\partial\theta}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta.
</math></center>

Integrando la primera ecuación respecto de <math>\rho</math> y comprobando la compatibilidad con la segunda, llegamos a

<center><math>
\psi(\rho,\theta)
= \left(\rho-\frac{1}{\rho}\right)\sin\theta
-\frac{1}{4\pi}\ln\rho + C,
</math></center>

donde <math>C</math> es una constante. Las curvas <math>\psi=\text{cte}</math> representan las líneas de corriente del nuevo flujo. Su representación numérica muestra cómo la introducción del término angular modifica el patrón de las corrientes, introduciendo una ligera asimetría respecto al caso inicial.

[[Archivo:lineasdecorrienteparaelnuevopotencial.png|thumb|400px|right|Líneas de corriente para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.3: Líneas de corriente del nuevo flujo

r = linspace(1,5,80);
a = linspace(0,2*pi,80);
[R,A] = meshgrid(r,a);

X = R.*cos(A);
Y = R.*sin(A);

% Campo de velocidades (reutilizamos las expresiones de 9.1)
u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

% Función de corriente psi
psi = (R - 1./R).*sin(A) - (1/(4*pi))*log(R);

figure
contourf(X,Y,psi,30,'LineStyle','-');
hold on
colorbar
quiver(X(1:3:end,1:3:end),Y(1:3:end,1:3:end), ...
Cx(1:3:end,1:3:end),Cy(1:3:end,1:3:end),1.3,'w');
plot(cos(a),sin(a),'k','LineWidth',2)
axis equal
axis([-4 4 -4 4])
title('Líneas de corriente para el nuevo potencial')
xlabel('EJE X'); ylabel('EJE Y');
hold off
}}

==Poster==
En el poster a continuación se muestra resumidamente el trabajo.
<center>[[Archivo:posterg68.JPG]]</center>

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-07T23:02:48Z

Jeremy.garcia.herrera: /* Representación de la función potencial y del campo de velocidades */

{{ TrabajoED | Flujo alrededor de un obstáculo circular. Grupo 68| [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] | Nicolás Rus Pérez, Jeremy García Herrera, Antonio Rodriguez Montes, Franco Luis Santillan, Juan Jose Salhua Palma.}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]
Un fluido es un medio contínuo deformable con la propiedad de que no transmite tensiones tangenciales. Estos no tienen una forma definida y por tanto adoptan la forma del recipiente en el que se encuentran, y pueden encontrarse tanto en estado líquido (agua) como gaseoso (aire).

En este trabajo se estudiará el flujo de un fluido alrededor de un obstáculo circular considerando a este un fluido incompresible. Esta propiedad indica que el fluido no disminuye su volumen al serle aplicado una fuerza, y es propia de los fluidos en estado líquido.

Por conveniencia se trabajará en coordenadas con base cilíndrica (polares).

==Mallado==
Se genera a continuación un mallado en el que se representan los puntos interiores a la región ocupada por el fluido. Por convenio se representará el obstáculo circular centrado en el origen de coordenadas y con radio 1. De esta forma la región ocupada por el fluido será todo el exterior de la circunferencia unidad.

El objetivo de realizar este mallado es dividir el medio continuo (la región ocupada por el fluido) en subdominios más pequeños denominados elementos finitos. El conjunto de estos elementos forman el mallado que cubre toda la región y permite realizar cálculos y simulaciones numéricas particulares en cada uno de ellos.

Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo <math>[-4,4]×[-4,4]</math>.
[[Archivo:Mallado1g68.JPG|thumb|right|Región del fluido y obstáculo]]
{{matlab|codigo=

r=linspace(1,5,50); %Definen el rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Genera el mallado

hold on
X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);
Z=0.*R;
mesh(X,Y,0*Z); %Dibuja la región
plot(cos(a),sin(a),'k','lineWidth',1); %Dibuja el obstáculo
axis([-4 4,-4 4]); %Se fijan los ejes
title ('Regíon Fluido'); %Definen los ejes y el nombre de la gráfica
xlabel 'EJE X'
ylabel 'EJE Y'
hold off
}}

==Campo de velocidad de las partículas==
Una vez definido el mallado de trabajo del fluido y el obstáculo, se procede a calcular las velocidades de las partículas del fluido en los distintos puntos de la región ocupada por el fluido incompresible. Partimos sabiendo que el campo de velocidades viene definido por el gradiente de la siguiente función potencial:

<center><math>\varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) </math></center>

A continuación se representa graficamente la función potencial:
[[Archivo:Mallado222g68.JPG|thumb|400px|right|Función potencial en la región del fluido]]
{{matlab|codigo=

r=linspace(1,5,30); %Define rho (r) y theta (a)
a=linspace(0,2*pi,30);

[R,A]=meshgrid(r,a); %CGenra el mallado

hold on
X=R.*cos(A); %Parametrización la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Define la función potencial
Z=f(R,A);
surf(X,Y,Z); %Representa la función
plot(1*cos(a),1*sin(a),'k','lineWidth',3); %Representa el obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de colores, define el titulo de la gráfica y los ejes
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
hold off
}}
Para calcular el campo de velocidades de las partículas del fluido es necesario el campo gradiente de la funcion potencial anterior, y dicho gradiente viene definido por la siguiente expresión:

<center>\(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)</center>

Aplicando dicha expresión a la función potencial, el gradiente resultante sera:

<center><math> \vec u(\vec e_\rho,\vec e_\theta,\vec e_z)=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta)]\vec e_\theta </math></center>

Dado que la representación en MATLAB se realiza en coordenadas cartesianas <math>(\vec i,\vec j,\vec k)</math>, y el gradiente calculado se encuentra en coordenadas ciliíndricas <math>(\vec e_\rho,\vec e_\theta,\vec e_z)</math> ,debemos de pasar el gradiente de la función potencial a las primeras. Esto es posible gracias a la matriz de cambio de base entre coordenadas cilíndricas y cartesianas, por lo que la multiplicamos por el vector gradiente dado <math>(\vec u)</math> en coordenadas cilíndricas.

Matriz de cambio de base de coordenadas cilindricas a coordenadas cartesianas:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z\\
\end{pmatrix}
</center>

Por lo que:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
(1-\frac{1}{\rho^2})\cos(\theta)\\
-(1+\frac{1}{\rho^2})\sin(\theta)\\
0\\
\end{pmatrix}
</center>
Resolviendo a coordenadas cartesianas:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta)\right) \vec j </math></center>

A continuación representamos el campo de velocidades de las partículas sobre las superficies equipotenciales del flujo del fluido.
[[Archivo:Mallado3g68.JPG|thumb|500px|right|Campo de velocidades de las partículas sobre las líneas equipotenciales del fluido]]
{{matlab|codigo=

r=linspace(1,5,50); %Define rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Crea la malla

X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Definimos la función potencial
Z=f(R,A); %Aplica la función
contour(X,Y,Z,30); %Dibuja las curvas de nivel
hold on %Define las componentes X e Y del gradiente
Cx=(1-(1./R.^2)).*cos(A).^2+(1+(1./R.^2)).*sin(A).^2;
Cy=(1-(1./R.^2)).*sin(A).*cos(A)-(1+(1./R.^2)).*sin(A).*cos(A);
quiver(X,Y,Cx,Cy); %Dibuja el campo de velocidades
plot(1*cos(a),1*sin(a),'k','lineWidth',1); %Representa del obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de color, titulo y nombre a los ejes
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}
Gracias a la gráfica anterior se observa como los vectores de velocidad definidos por <math> \vec u </math> son siempre perpendiculares a las líneas equipotenciales del fluido.
<center>[[Archivo:Mallado33.JPG]]</center>

==Rotacional Nulo y Divergencia Nula==
===Rotacional Nulo===
El rotacional nulo del campo de velocidades indica que el flujo es '''irrotacional''' , es decir , carece de vórtices o giros locales. Esto significa que una partícula de fluido no rota sobre sí misma mientras se desplaza, característica típica de los flujos potenciales.

El operador rotacional nos indica la dirección y la velocidad del giro del campo en cada punto.

A continuación se demuestra que el rotacional es nulo:

<center><math>
\nabla\times \vec{u} = \frac{1}{\rho} \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & \rho(-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix}=</math>

<math>
= \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} \left[ \rho \left( -(1 + \frac{1}{{\rho^2}})\sin(\theta) \right) \right]\vec{e}_z - \frac{\partial}{\partial \theta} \left[ \cos(\theta) \left(1 - \frac{1}{\rho^2}\right) \right] \vec{e}_z \right]=\frac{1}{\rho} \left[ -(1 - \frac{1}{\rho^2})\sin(\theta) + \sin(\theta)(1 - \frac{1}{\rho^2}) \right]=0
</math>

<math>\boxed { \nabla \times \bar { u } =0 } </math></center>

Se obtiene que el rotacional es 0, lo cual indica que las partículas del fluido no están girando

===Divergencia Nula===
Interpretación de la divergencia nula: La condición ∇·u = 0 corresponde a un flujo incompresible, donde la densidad del fluido permanece constante y el volumen de cualquier elemento material no cambia a lo largo del tiempo. Esto refleja la conservación local de la masa en ausencia de fuentes o sumideros.

<center><math>
\nabla \cdot \vec{u}(\rho, \theta, z) = \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} (\rho u_\rho) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (\rho u_z) \right\}= </math>

<math>
= \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} \left( (\rho - \frac{1}{\rho})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) \right) \right\}=
</math>
<math>
\frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0
</math>

<math>\boxed { div(\bar { u } )=0 } </math></center>

== Líneas de corriente ==
Las líneas de corriente del flujo alrededor de un cilindro circular (en sección transversal) presentan un patrón simétrico con respecto al eje horizontal (eje x) y antisimétrico con respecto al eje vertical (eje y).
[[Archivo:Linea_corrienteTangencial.png|thumb|400px|right|Líneas de corrientes]]

Las líneas de corriente <math>\psi = \text{cte}</math> describen las trayectorias del flujo potencial alrededor del obstáculo circular, mostrando simetría y adaptación suave a la geometría. Estas líneas son tangentes al campo de velocidades en cada punto, verificando que <math>\vec{u}\cdot\nabla\psi = 0</math> . El patrón obtenido revela aceleración del fluido cerca del cilindro, condición de no penetración en la superficie, y recuperación del flujo uniforme lejos del obstáculo, características típicas de un flujo incompresible e irrotacional.

<math>\vec{v}</math> se calculará con : <math>\vec{v}=\vec{k}\times \vec{u}</math>

<math>\vec{v}</math> es irrotacional dadon que el rotacional calculado en el apartado 3.1 es nulo.
De igual forma iene un potencial escalar <math>\psi </math> el cual tiene <math>\vec{v}</math> de gradiente.

Para el potencial escalar <math>\psi </math> se resolverá este sistema de ecuaciones: <math>\nabla\psi =\vec{v}</math> cumpliendo la condición anterior de líneas <math>\psi =cte</math>.

Empezamos las operaciones:
<math>\vec{v} =\vec{k}\times \vec{u} = \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
0 & 0 & 1 \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & (-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

<center><math>\boxed{\frac{\partial \psi}{\partial \rho} = V\rho,
\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = V\theta,
\frac{\partial \psi}{\partial z} = Vz}
</math></center>

<math>\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho = \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

<math>\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})\Rightarrow\psi(\rho,\theta,Z)= \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

Dando como resultado
<math>\psi: \sin{\theta}(\rho-\frac{1}{\rho})=cte</math>
[[Archivo:Psi_lineasCorriente.png|thumb|400px|right|Líneas de corrientes]]
{{matlab|codigo=
% Definir dominio en coordenadas cilíndricas
r = linspace(1, 5, 50); % Coordenada radial (ρ) desde r=1 (obstáculo)
a = linspace(0, 2*pi, 50); % Coordenada angular (θ)

% Crear malla
[R, A] = meshgrid(r, a);

% Convertir a coordenadas cartesianas
X = R .* cos(A);
Y = R .* sin(A);

% 1. FUNCIÓN POTENCIAL φ (phi)
phi = (R + (1./R)) .* cos(A);
% 2. CAMPO DE VELOCIDADES u = ∇φ

% Componentes en coordenadas cilíndricas:
% u_ρ = (1 - 1/ρ²)cos(θ)
% u_θ = -(1 + 1/ρ²)sin(θ)

u_rho = (1 - (1./R.^2)) .* cos(A);
u_theta = -(1 + (1./R.^2)) .* sin(A);

% Convertir a coordenadas cartesianas
u_x = u_rho .* cos(A) - u_theta .* sin(A);
u_y = u_rho .* sin(A) + u_theta .* cos(A);

% 3. FUNCIÓN DE CORRIENTE ψ (psi)
psi = (R - 1./R) .* sin(A);

% 4. GRÁFICAS

figure('Position', [100, 100, 1200, 800]);

% Graf 1: Líneas de corriente (ψ) y campo de velocidades
subplot(2, 2, 2);
contourf(X, Y, psi, 30, 'LineStyle', '-', 'LineWidth', 0.5);
hold on;
colorbar;
title('Líneas de Corriente (ψ) y Campo de Velocidades');
xlabel('EJE X');
ylabel('EJE Y');
axis equal;
axis([-4, 4, -4, 4]);
grid on;

% Dibujar campo de velocidades
quiver(X(1:3:end, 1:3:end), Y(1:3:end, 1:3:end), ...
u_x(1:3:end, 1:3:end), u_y(1:3:end, 1:3:end), ...
1.5, 'w', 'LineWidth', 0.7, 'MaxHeadSize', 0.5);

% Dibujar obstáculo circular
plot(cos(theta_circ), sin(theta_circ), 'k', 'LineWidth', 2, ...
'Color', [0.2, 0.2, 0.2]);
fill(cos(theta_circ), sin(theta_circ), [0.8, 0.8, 0.8], ...
'EdgeColor', [0.2, 0.2, 0.2], 'LineWidth', 1.5);
hold off;
}}

== Velocidad extrema y puntos de remanso==

En la frontera del obstáculo, definida por la circunferencia unidad <math>\rho=1</math>, el módulo de la velocidad viene dado por

<center><math>|\vec u| = 2|\sin\theta|</math></center>

De esta expresión se deduce que la velocidad máxima se alcanza cuando

<math>\sin\theta=\pm1</math>, es decir, en

<center><math>\theta=\frac{\pi}{2},\ \frac{3\pi}{2}</math></center>

que corresponden a los puntos cartesianos

<center><math>(0,1)\ \text{y}\ (0,-1)</math></center>

Estos son los puntos donde la velocidad del fluido es máxima sobre la frontera.

Por otro lado, la velocidad se anula cuando

<center><math>\sin\theta=0 \Rightarrow \theta=0,\ \pi</math></center>

lo que corresponde a los puntos

<center><math>(1,0)\ \text{y}\ (-1,0)</math></center>

Estos puntos reciben el nombre de '''puntos de remanso'''.

Para verificarlo numéricamente en MATLAB se ha representado el módulo de la velocidad sobre la frontera:
[[Archivo:Subir1.png|thumb|500px|center|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular. La velocidad es nula en los puntos de remanso y máxima en la dirección vertical.]]
{{matlab|codigo=
nr=100; nt=200;
r=linspace(1,5,nr);
th=linspace(0,2*pi,nt);
[R,TH]=meshgrid(r,th);

ur=(1-1./R.^2).*cos(TH);
uth=-(1+1./R.^2).*sin(TH);

Ux=ur.*cos(TH)-uth.*sin(TH);
Uy=ur.*sin(TH)+uth.*cos(TH);

Umag=sqrt(Ux.^2+Uy.^2);

UmagS=Umag(:,1);
THS=TH(:,1);

[Umax,imax]=max(UmagS);
[Umin,imin]=min(UmagS);

thetaMax=THS(imax);
thetaMin=THS(imin);

figure
plot(THS,UmagS), grid on
hold on
plot(thetaMax,Umax,'ro',thetaMin,Umin,'bo')
title('Modulo de la velocidad sobre la frontera S')
xlabel('\theta'), ylabel('|\itu|')
hold off
}}

== Presión del fluido y ecuación de Bernoulli==

Se supone que la densidad del fluido es constante e igual a

<center><math>\rho=2</math></center>

y se cumple la ecuación de Bernoulli

<center><math>\frac12\rho|\vec u|^2+p=C</math></center>

Despejando la presión se obtiene

<center><math>p=C-\frac12\rho|\vec u|^2</math></center>

Tomando para la constante el valor <math>C=1</math>, la presión queda completamente determinada en toda la región del fluido.

El campo de presiones se ha representado en MATLAB mediante:

{{matlab|codigo=
rho=2;
C=1;

P=C-0.5*rho*Umag.^2;

PS=P(:,1);
[Pmax,ipmax]=max(PS);
[Pmin,ipmin]=min(PS);

thetaPmax=THS(ipmax);
thetaPmin=THS(ipmin);

figure
contourf(R.*cos(TH),R.*sin(TH),P,30,'LineStyle','none');
colorbar
hold on
viscircles([0 0],1,'Color','k')
axis equal
title('Presión p(x,y)')
xlabel('x'), ylabel('y')
hold off
}}
{| class="wikitable" style="border:0; margin:auto; background:none;"
| style="border:0;" | [[Archivo:Subir3.png|thumb|500px|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular]]
| style="border:0;" | [[Archivo:Subir2.png|thumb|500px|Módulo de la velocidad ∣𝑢⃗∣ sobre la frontera del obstáculo 𝑟=1 en función del ángulo 𝜃]]
|}
Se observa que la presión es máxima en los puntos de remanso, donde la velocidad es nula, y mínima en los puntos donde la velocidad alcanza su valor máximo. Esto concuerda con la ecuación de Bernoulli, ya que un aumento del módulo de la velocidad implica necesariamente una disminución de la presión.

==Velocidad y presión al rodear el obstáculo.==

Ya teniendo el campo de presiones se puede apreciar que tiene relación con la velocidad, ya que cuando disminuye la presión aumenta el módulo de la velocidad, tratándose de un resultado lógico debido a que cuando hablamos de un fluido incompresible su valor de Trinomio de Bernoulli ha de ser constante.
Lo podemos ver mediante la siguiente expresión: (Ec. Bernoulli).

<math>B=z+P/γ+v^2/2g</math>

B=valor del trinomio de Bernoulli

P= presión

z= altura

v=velocidad que lleva el fluido en determinado punto

g=valor de la gravedad

γ=valor del peso específico

Donde podemos ver que si la presión aumenta la velocidad tiene que disminuir y viceversa.

==Paradoja de D'Alembert.==

La magnitud pn representa la fuerza que el fluido ejerce de forma puntual sobre la superficie del obstáculo. Para estudiar si existe empuje neto en la dirección horizontal, basta analizar la suma de las proyecciones de todas esas fuerzas sobre el vector i. Si dicha suma resulta nula, concluiremos que no hay fuerza total en esa dirección.

Del apartado anterior conocemos ya la presión en la frontera del círculo unidad:
<math>p=10-∥u∥^2=9-4 sin⁡^2 θ-4sin⁡θ</math>

Cálculo del vector normal

Sobre la circunferencia de radio 1 el vector normal exterior es

<math>n=-eρ.</math>

Para expresarlo en coordenadas cartesianas utilizamos la matriz de cambio de base de polares a cartesianas. Tomando su traspuesta se obtiene:

\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}

Aplicando esta matriz a <math>(1,0,0)^T</math>, que corresponde al versor radial eρ, se obtiene:
\begin{pmatrix}
\cos \theta &\\
sen \theta &\\
0
\end{pmatrix}
<math>=cos⁡θi+sin⁡θj</math>

Al tomar el signo negativo para el normal exterior,

<math>n=cos⁡θi-sin⁡θj</math>,

de donde su proyección sobre i es:

<math>n⋅i=-cos⁡θ</math>.

Integración sobre la frontera

El borde del dominio se parametriza como
<math>(ρ,θ)=(1,θ), con θ∈[0,2π)</math>.

La contribución total de la presión en la dirección horizontal es:
<math>\int_{0}^{2π}p(n⋅i)dθ=\int_{0}^{2π}p(-cos⁡θ)dθ.</math>

Sustituyendo la expresión de p:
<math>\int_{0}^{2π}(9-4sin⁡^2 θ-4sin⁡θ)(-cos⁡θ)dθ.</math>

Integramos término a término:
<math>\int_{0}^{2π}-9cos⁡θdθ+\int_{0}^{2π}4sin⁡^2 θcos⁡θdθ+\int_{0}^{2π}4sin⁡θcos⁡θdθ.</math>

Cada una de estas integrales es periódica con una primitiva continua en todo el intervalo; evaluándolas entre 0 y 2π se obtiene:
<math>[-9sin⁡θ]_{0}^{2π}=0,[(4sin⁡^3 θ)/3]_{0}^{2π}=0,[2sin⁡^2 θ]_{0}^{2π}=0.</math>

Por tanto,
<math>\int_{0}^{2π}p(n⋅i)dθ=0.</math>

Conclusión

La componente horizontal de la fuerza total ejercida por el fluido sobre el obstáculo es nula. En consecuencia, el fluido no produce empuje neto en la dirección i.

==Cambio de función potencial==

En este apartado repetimos los procesos de los puntos (2), (3) y (4), pero tomando ahora como función potencial

<center><math>\varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}</math></center>

El obstáculo sigue siendo la circunferencia unidad centrada en el origen, y el fluido ocupa la región exterior, trabajándose de nuevo en coordenadas cilíndricas (polares).

===Representación de la función potencial y del campo de velocidades===

En primer lugar representamos la nueva función potencial <math>\varphi</math> sobre la región ocupada por el fluido, manteniendo el mismo mallado del apartado inicial. Esta representación permite visualizar cómo el término angular <math>\theta/(4\pi)</math> introduce una ligera inclinación en la superficie potencial respecto al caso anterior.

[[Archivo:nuevafuncionpotencial.png|thumb|400px|right|Nueva función potencial del flujo con término angular]]

La velocidad de las partículas del fluido viene dada de nuevo por el gradiente del potencial. En coordenadas cilíndricas, el campo de velocidades es

<center><math>\vec u(\rho,\theta,z)=\nabla\varphi
= u_\rho\,\vec e_\rho+u_\theta\,\vec e_\theta+u_z\,\vec e_z</math></center>

donde

<center><math>
u_\rho=\frac{\partial\varphi}{\partial\rho}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta, \qquad
u_\theta=\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}
= -\left(1+\frac{1}{\rho^2}\right)\sin\theta+\frac{1}{4\pi\rho},\qquad
u_z=0.
</math></center>

Para poder representarlo en MATLAB se pasa el campo a coordenadas cartesianas mediante la matriz de cambio de base de cilíndricas a cartesianas:

<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k
\end{pmatrix}
=
\begin{pmatrix}
\cos\theta & -\sin\theta & 0\\
\sin\theta & \cos\theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z
\end{pmatrix}
</center>

Aplicando esta matriz al vector <math>(u_\rho,u_\theta,0)^T</math> se obtiene

<center><math>
\vec u(\vec i,\vec j,\vec k)
=\left[(1-\tfrac{1}{\rho^2})\cos^2\theta+(1+\tfrac{1}{\rho^2})\sin^2\theta
-\frac{\sin\theta}{4\pi\rho}\right]\vec i
+\left[(1-\tfrac{1}{\rho^2})\sin\theta\cos\theta
-(1+\tfrac{1}{\rho^2})\sin\theta\cos\theta
+\frac{\cos\theta}{4\pi\rho}\right]\vec j.
</math></center>

Sobre las curvas de nivel de <math>\varphi</math> se representa este nuevo campo de velocidades:

[[Archivo:campodevelocidadesnuevopotencial.png|thumb|500px|right|Campo de velocidades sobre las líneas equipotenciales para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.1: Potencial y campo de velocidades para la nueva función φ

r = linspace(1,5,30); % Coordenada radial (rho)
a = linspace(0,2*pi,30); % Coordenada angular (theta)
[R,A] = meshgrid(r,a); % Mallado en (rho,theta)

% Conversión a coordenadas cartesianas
X = R.*cos(A);
Y = R.*sin(A);

% Nueva función potencial
phi = (R + 1./R).*cos(A) + A/(4*pi);

figure
surf(X,Y,phi);
shading interp
hold on
plot(1*cos(a),1*sin(a),'k','LineWidth',2); % Obstáculo
axis([-4 4 -4 4])
title('Nueva función potencial')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off

% Campo de velocidades (en cartesianas) sobre curvas de nivel de φ
figure
contour(X,Y,phi,30);
hold on

u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

quiver(X,Y,Cx,Cy)
plot(1*cos(a),1*sin(a),'k','LineWidth',1)
axis equal
axis([-4 4 -4 4])
title('Campo de velocidades (nuevo potencial)')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off
}}

Gracias a la representación anterior se observa que los vectores velocidad siguen siendo ortogonales a las líneas equipotenciales del nuevo potencial.

===Rotacional nulo y divergencia nula del nuevo campo===

Al ser <math>\vec u</math> el gradiente de un potencial escalar, su rotacional sigue siendo nulo:

<center><math>\nabla\times\vec u = \vec 0.</math></center>

En coordenadas cilíndricas se verifica explícitamente mediante la expresión

<center><math>
\nabla\times\vec u =
\frac{1}{\rho}
\begin{vmatrix}
\vec e_\rho & \vec e_\theta & \vec e_z\\
\frac{\partial}{\partial\rho} &
\frac{\partial}{\partial\theta} &
\frac{\partial}{\partial z}\\
u_\rho & \rho u_\theta & 0
\end{vmatrix}
= \vec 0,
</math></center>

lo que muestra que el flujo sigue siendo irrotacional.

Por otro lado, la divergencia del nuevo campo de velocidades viene dada por

<center><math>
\nabla\cdot\vec u(\rho,\theta,z)
=\frac{1}{\rho}\left\{
\frac{\partial}{\partial\rho}(\rho u_\rho)
+\frac{\partial}{\partial\theta}(u_\theta)
+\frac{\partial}{\partial z}(\rho u_z)
\right\},
</math></center>

y sustituyendo las expresiones de <math>u_\rho</math> y <math>u_\theta</math> se obtiene de nuevo

<center><math>\nabla\cdot\vec u = 0.</math></center>

La divergencia nula indica que el flujo sigue siendo incompresible: no hay creación ni destrucción de volumen de fluido en el dominio.

===Líneas de corriente del nuevo flujo===

Para determinar las líneas de corriente asociadas al nuevo campo de velocidades, consideramos el campo <math>\vec v</math> ortogonal a <math>\vec u</math>, definido por

<center><math>\vec v = \vec k \times \vec u.</math></center>

En coordenadas cilíndricas se tiene

<center>
<math>
\vec v
= \left(\left(1+\tfrac{1}{\rho^2}\right)\sin\theta - \tfrac{1}{4\pi\rho}\right)\vec e_\rho
+ \left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta.
</math>
</center>

Dado que el rotacional de <math>\vec v</math> también es nulo, existe un potencial escalar <math>\psi</math> tal que

<center><math>\nabla\psi = \vec v.</math></center>

A partir de las componentes de <math>\vec v</math> se obtiene el siguiente sistema:

<center><math>
\frac{\partial\psi}{\partial\rho}
= \left(1+\frac{1}{\rho^2}\right)\sin\theta - \frac{1}{4\pi\rho},\qquad
\frac{1}{\rho}\frac{\partial\psi}{\partial\theta}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta.
</math></center>

Integrando la primera ecuación respecto de <math>\rho</math> y comprobando la compatibilidad con la segunda, llegamos a

<center><math>
\psi(\rho,\theta)
= \left(\rho-\frac{1}{\rho}\right)\sin\theta
-\frac{1}{4\pi}\ln\rho + C,
</math></center>

donde <math>C</math> es una constante. Las curvas <math>\psi=\text{cte}</math> representan las líneas de corriente del nuevo flujo. Su representación numérica muestra cómo la introducción del término angular modifica el patrón de las corrientes, introduciendo una ligera asimetría respecto al caso inicial.

[[Archivo:LineasCorriente9_g68.JPG|thumb|400px|right|Líneas de corriente para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.3: Líneas de corriente del nuevo flujo

r = linspace(1,5,80);
a = linspace(0,2*pi,80);
[R,A] = meshgrid(r,a);

X = R.*cos(A);
Y = R.*sin(A);

% Campo de velocidades (reutilizamos las expresiones de 9.1)
u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

% Función de corriente psi
psi = (R - 1./R).*sin(A) - (1/(4*pi))*log(R);

figure
contourf(X,Y,psi,30,'LineStyle','-');
hold on
colorbar
quiver(X(1:3:end,1:3:end),Y(1:3:end,1:3:end), ...
Cx(1:3:end,1:3:end),Cy(1:3:end,1:3:end),1.3,'w');
plot(cos(a),sin(a),'k','LineWidth',2)
axis equal
axis([-4 4 -4 4])
title('Líneas de corriente para el nuevo potencial')
xlabel('EJE X'); ylabel('EJE Y');
hold off
}}

==Poster==
En el poster a continuación se muestra resumidamente el trabajo.
<center>[[Archivo:posterg68.JPG]]</center>

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-07T22:58:34Z

Jeremy.garcia.herrera:

{{ TrabajoED | Flujo alrededor de un obstáculo circular. Grupo 68| [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] | Nicolás Rus Pérez, Jeremy García Herrera, Antonio Rodriguez Montes, Franco Luis Santillan, Juan Jose Salhua Palma.}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]
Un fluido es un medio contínuo deformable con la propiedad de que no transmite tensiones tangenciales. Estos no tienen una forma definida y por tanto adoptan la forma del recipiente en el que se encuentran, y pueden encontrarse tanto en estado líquido (agua) como gaseoso (aire).

En este trabajo se estudiará el flujo de un fluido alrededor de un obstáculo circular considerando a este un fluido incompresible. Esta propiedad indica que el fluido no disminuye su volumen al serle aplicado una fuerza, y es propia de los fluidos en estado líquido.

Por conveniencia se trabajará en coordenadas con base cilíndrica (polares).

==Mallado==
Se genera a continuación un mallado en el que se representan los puntos interiores a la región ocupada por el fluido. Por convenio se representará el obstáculo circular centrado en el origen de coordenadas y con radio 1. De esta forma la región ocupada por el fluido será todo el exterior de la circunferencia unidad.

El objetivo de realizar este mallado es dividir el medio continuo (la región ocupada por el fluido) en subdominios más pequeños denominados elementos finitos. El conjunto de estos elementos forman el mallado que cubre toda la región y permite realizar cálculos y simulaciones numéricas particulares en cada uno de ellos.

Para ilustrar que el fluido ocupa el exterior de un circulo, se dibujarán los ejes en el intervalo <math>[-4,4]×[-4,4]</math>.
[[Archivo:Mallado1g68.JPG|thumb|right|Región del fluido y obstáculo]]
{{matlab|codigo=

r=linspace(1,5,50); %Definen el rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Genera el mallado

hold on
X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);
Z=0.*R;
mesh(X,Y,0*Z); %Dibuja la región
plot(cos(a),sin(a),'k','lineWidth',1); %Dibuja el obstáculo
axis([-4 4,-4 4]); %Se fijan los ejes
title ('Regíon Fluido'); %Definen los ejes y el nombre de la gráfica
xlabel 'EJE X'
ylabel 'EJE Y'
hold off
}}

==Campo de velocidad de las partículas==
Una vez definido el mallado de trabajo del fluido y el obstáculo, se procede a calcular las velocidades de las partículas del fluido en los distintos puntos de la región ocupada por el fluido incompresible. Partimos sabiendo que el campo de velocidades viene definido por el gradiente de la siguiente función potencial:

<center><math>\varphi (\rho ,\theta)=(\rho +\frac{1}{\rho})\cos (\theta ) </math></center>

A continuación se representa graficamente la función potencial:
[[Archivo:Mallado222g68.JPG|thumb|400px|right|Función potencial en la región del fluido]]
{{matlab|codigo=

r=linspace(1,5,30); %Define rho (r) y theta (a)
a=linspace(0,2*pi,30);

[R,A]=meshgrid(r,a); %CGenra el mallado

hold on
X=R.*cos(A); %Parametrización la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Define la función potencial
Z=f(R,A);
surf(X,Y,Z); %Representa la función
plot(1*cos(a),1*sin(a),'k','lineWidth',3); %Representa el obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de colores, define el titulo de la gráfica y los ejes
title ('Función potencial');
xlabel ('EJE X');
ylabel ('EJE Y');
hold off
}}
Para calcular el campo de velocidades de las partículas del fluido es necesario el campo gradiente de la funcion potencial anterior, y dicho gradiente viene definido por la siguiente expresión:

<center>\(\nabla f(\rho, \theta, z) = \frac{\partial f}{\partial \rho} \vec e_\rho + \frac{1}{\rho} \frac{\partial f}{\partial \theta} \vec e_\theta + \frac{\partial f}{\partial z} \vec e_z\)</center>

Aplicando dicha expresión a la función potencial, el gradiente resultante sera:

<center><math> \vec u(\vec e_\rho,\vec e_\theta,\vec e_z)=\nabla \varphi=[(1-\frac{1}{\rho^2})\cos(\theta)]\vec e_\rho-[(1+\frac{1}{\rho^2})\sin(\theta)]\vec e_\theta </math></center>

Dado que la representación en MATLAB se realiza en coordenadas cartesianas <math>(\vec i,\vec j,\vec k)</math>, y el gradiente calculado se encuentra en coordenadas ciliíndricas <math>(\vec e_\rho,\vec e_\theta,\vec e_z)</math> ,debemos de pasar el gradiente de la función potencial a las primeras. Esto es posible gracias a la matriz de cambio de base entre coordenadas cilíndricas y cartesianas, por lo que la multiplicamos por el vector gradiente dado <math>(\vec u)</math> en coordenadas cilíndricas.

Matriz de cambio de base de coordenadas cilindricas a coordenadas cartesianas:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z\\
\end{pmatrix}
</center>

Por lo que:
<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k\\
\end{pmatrix}
=
\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
(1-\frac{1}{\rho^2})\cos(\theta)\\
-(1+\frac{1}{\rho^2})\sin(\theta)\\
0\\
\end{pmatrix}
</center>
Resolviendo a coordenadas cartesianas:
<center><math> \vec u (\vec i,\vec j,\vec k) =\nabla \varphi=\left( (1 - \frac{1}{{\rho^2}}) \cdot \cos^2(\theta) + (1 + \frac{1}{{\rho^2}}) \cdot \sin^2(\theta)\right) \vec i + \left( (1 - \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta) - (1 + \frac{1}{{\rho^2}}) \cdot \sin(\theta) \cdot \cos(\theta)\right) \vec j </math></center>

A continuación representamos el campo de velocidades de las partículas sobre las superficies equipotenciales del flujo del fluido.
[[Archivo:Mallado3g68.JPG|thumb|500px|right|Campo de velocidades de las partículas sobre las líneas equipotenciales del fluido]]
{{matlab|codigo=

r=linspace(1,5,50); %Define rho (r) y theta (a)
a=linspace(0,2*pi,50);

[R,A]=meshgrid(r,a); %Crea la malla

X=R.*cos(A); %Parametriza la superficie
Y=R.*sin(A);

f=@(r,a)(r+(1./r)).*cos(a); %Definimos la función potencial
Z=f(R,A); %Aplica la función
contour(X,Y,Z,30); %Dibuja las curvas de nivel
hold on %Define las componentes X e Y del gradiente
Cx=(1-(1./R.^2)).*cos(A).^2+(1+(1./R.^2)).*sin(A).^2;
Cy=(1-(1./R.^2)).*sin(A).*cos(A)-(1+(1./R.^2)).*sin(A).*cos(A);
quiver(X,Y,Cx,Cy); %Dibuja el campo de velocidades
plot(1*cos(a),1*sin(a),'k','lineWidth',1); %Representa del obstáculo
axis([-4,4,-4,4]);
colorbar; %Añade la barra de color, titulo y nombre a los ejes
title ('Campo de velocidades');
xlabel ('EJE X');
ylabel ('EJE Y');
axis equal
hold off
}}
Gracias a la gráfica anterior se observa como los vectores de velocidad definidos por <math> \vec u </math> son siempre perpendiculares a las líneas equipotenciales del fluido.
<center>[[Archivo:Mallado33.JPG]]</center>

==Rotacional Nulo y Divergencia Nula==
===Rotacional Nulo===
El rotacional nulo del campo de velocidades indica que el flujo es '''irrotacional''' , es decir , carece de vórtices o giros locales. Esto significa que una partícula de fluido no rota sobre sí misma mientras se desplaza, característica típica de los flujos potenciales.

El operador rotacional nos indica la dirección y la velocidad del giro del campo en cada punto.

A continuación se demuestra que el rotacional es nulo:

<center><math>
\nabla\times \vec{u} = \frac{1}{\rho} \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
\frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & \rho(-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix}=</math>

<math>
= \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} \left[ \rho \left( -(1 + \frac{1}{{\rho^2}})\sin(\theta) \right) \right]\vec{e}_z - \frac{\partial}{\partial \theta} \left[ \cos(\theta) \left(1 - \frac{1}{\rho^2}\right) \right] \vec{e}_z \right]=\frac{1}{\rho} \left[ -(1 - \frac{1}{\rho^2})\sin(\theta) + \sin(\theta)(1 - \frac{1}{\rho^2}) \right]=0
</math>

<math>\boxed { \nabla \times \bar { u } =0 } </math></center>

Se obtiene que el rotacional es 0, lo cual indica que las partículas del fluido no están girando

===Divergencia Nula===
Interpretación de la divergencia nula: La condición ∇·u = 0 corresponde a un flujo incompresible, donde la densidad del fluido permanece constante y el volumen de cualquier elemento material no cambia a lo largo del tiempo. Esto refleja la conservación local de la masa en ausencia de fuentes o sumideros.

<center><math>
\nabla \cdot \vec{u}(\rho, \theta, z) = \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} (\rho u_\rho) + \frac{\partial}{\partial \theta} (u_\theta) + \frac{\partial}{\partial z} (\rho u_z) \right\}= </math>

<math>
= \frac{1}{\rho} \left\{ \frac{\partial}{\partial \rho} \left( (\rho - \frac{1}{\rho})\cos(\theta) \right) + \frac{\partial}{\partial \theta} \left( -\left(1 + \frac{1}{\rho^2}\right)\sin(\theta) \right) \right\}=
</math>
<math>
\frac{1}{\rho} \left\{ (1 + \frac{1}{\rho^2})\cos(\theta) -\left(1 + \frac{1}{\rho^2}\right)\cos(\theta) \right\}= 0
</math>

<math>\boxed { div(\bar { u } )=0 } </math></center>

== Líneas de corriente ==
Las líneas de corriente del flujo alrededor de un cilindro circular (en sección transversal) presentan un patrón simétrico con respecto al eje horizontal (eje x) y antisimétrico con respecto al eje vertical (eje y).
[[Archivo:Linea_corrienteTangencial.png|thumb|400px|right|Líneas de corrientes]]

Las líneas de corriente <math>\psi = \text{cte}</math> describen las trayectorias del flujo potencial alrededor del obstáculo circular, mostrando simetría y adaptación suave a la geometría. Estas líneas son tangentes al campo de velocidades en cada punto, verificando que <math>\vec{u}\cdot\nabla\psi = 0</math> . El patrón obtenido revela aceleración del fluido cerca del cilindro, condición de no penetración en la superficie, y recuperación del flujo uniforme lejos del obstáculo, características típicas de un flujo incompresible e irrotacional.

<math>\vec{v}</math> se calculará con : <math>\vec{v}=\vec{k}\times \vec{u}</math>

<math>\vec{v}</math> es irrotacional dadon que el rotacional calculado en el apartado 3.1 es nulo.
De igual forma iene un potencial escalar <math>\psi </math> el cual tiene <math>\vec{v}</math> de gradiente.

Para el potencial escalar <math>\psi </math> se resolverá este sistema de ecuaciones: <math>\nabla\psi =\vec{v}</math> cumpliendo la condición anterior de líneas <math>\psi =cte</math>.

Empezamos las operaciones:
<math>\vec{v} =\vec{k}\times \vec{u} = \begin{vmatrix}
\vec{e}_\rho & \vec{e}_\theta & \vec{e}_z \\
0 & 0 & 1 \\
(1 - \frac{1}{{\rho^2}})\cos(\theta) & (-(1 + \frac{1}{{\rho^2}})\sin(\theta)) & 0
\end{vmatrix} = ((1+\frac{1}{\rho^2})\sin{\theta})\vec{e}_\rho + (1-\frac{1}{\rho^2})\cos{\theta}\vec{e}_\theta</math>

<center><math>\boxed{\frac{\partial \psi}{\partial \rho} = V\rho,
\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = V\theta,
\frac{\partial \psi}{\partial z} = Vz}
</math></center>

<math>\psi(\rho,\theta,Z) =\int((1+\frac{1}{\rho^2})\sin{\theta})\partial\rho = \sin{\theta}\int(1+\frac{1}{\rho^2})\partial\rho = \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

<math>\frac{1}{\rho}\frac{\partial \psi}{\partial \theta} = \cos{\theta}(\rho+\frac{1}{\rho})\Rightarrow\psi(\rho,\theta,Z)= \sin{\theta}(\rho-\frac{1}{\rho})+ C</math>

Dando como resultado
<math>\psi: \sin{\theta}(\rho-\frac{1}{\rho})=cte</math>
[[Archivo:Psi_lineasCorriente.png|thumb|400px|right|Líneas de corrientes]]
{{matlab|codigo=
% Definir dominio en coordenadas cilíndricas
r = linspace(1, 5, 50); % Coordenada radial (ρ) desde r=1 (obstáculo)
a = linspace(0, 2*pi, 50); % Coordenada angular (θ)

% Crear malla
[R, A] = meshgrid(r, a);

% Convertir a coordenadas cartesianas
X = R .* cos(A);
Y = R .* sin(A);

% 1. FUNCIÓN POTENCIAL φ (phi)
phi = (R + (1./R)) .* cos(A);
% 2. CAMPO DE VELOCIDADES u = ∇φ

% Componentes en coordenadas cilíndricas:
% u_ρ = (1 - 1/ρ²)cos(θ)
% u_θ = -(1 + 1/ρ²)sin(θ)

u_rho = (1 - (1./R.^2)) .* cos(A);
u_theta = -(1 + (1./R.^2)) .* sin(A);

% Convertir a coordenadas cartesianas
u_x = u_rho .* cos(A) - u_theta .* sin(A);
u_y = u_rho .* sin(A) + u_theta .* cos(A);

% 3. FUNCIÓN DE CORRIENTE ψ (psi)
psi = (R - 1./R) .* sin(A);

% 4. GRÁFICAS

figure('Position', [100, 100, 1200, 800]);

% Graf 1: Líneas de corriente (ψ) y campo de velocidades
subplot(2, 2, 2);
contourf(X, Y, psi, 30, 'LineStyle', '-', 'LineWidth', 0.5);
hold on;
colorbar;
title('Líneas de Corriente (ψ) y Campo de Velocidades');
xlabel('EJE X');
ylabel('EJE Y');
axis equal;
axis([-4, 4, -4, 4]);
grid on;

% Dibujar campo de velocidades
quiver(X(1:3:end, 1:3:end), Y(1:3:end, 1:3:end), ...
u_x(1:3:end, 1:3:end), u_y(1:3:end, 1:3:end), ...
1.5, 'w', 'LineWidth', 0.7, 'MaxHeadSize', 0.5);

% Dibujar obstáculo circular
plot(cos(theta_circ), sin(theta_circ), 'k', 'LineWidth', 2, ...
'Color', [0.2, 0.2, 0.2]);
fill(cos(theta_circ), sin(theta_circ), [0.8, 0.8, 0.8], ...
'EdgeColor', [0.2, 0.2, 0.2], 'LineWidth', 1.5);
hold off;
}}

== Velocidad extrema y puntos de remanso==

En la frontera del obstáculo, definida por la circunferencia unidad <math>\rho=1</math>, el módulo de la velocidad viene dado por

<center><math>|\vec u| = 2|\sin\theta|</math></center>

De esta expresión se deduce que la velocidad máxima se alcanza cuando

<math>\sin\theta=\pm1</math>, es decir, en

<center><math>\theta=\frac{\pi}{2},\ \frac{3\pi}{2}</math></center>

que corresponden a los puntos cartesianos

<center><math>(0,1)\ \text{y}\ (0,-1)</math></center>

Estos son los puntos donde la velocidad del fluido es máxima sobre la frontera.

Por otro lado, la velocidad se anula cuando

<center><math>\sin\theta=0 \Rightarrow \theta=0,\ \pi</math></center>

lo que corresponde a los puntos

<center><math>(1,0)\ \text{y}\ (-1,0)</math></center>

Estos puntos reciben el nombre de '''puntos de remanso'''.

Para verificarlo numéricamente en MATLAB se ha representado el módulo de la velocidad sobre la frontera:
[[Archivo:Subir1.png|thumb|500px|center|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular. La velocidad es nula en los puntos de remanso y máxima en la dirección vertical.]]
{{matlab|codigo=
nr=100; nt=200;
r=linspace(1,5,nr);
th=linspace(0,2*pi,nt);
[R,TH]=meshgrid(r,th);

ur=(1-1./R.^2).*cos(TH);
uth=-(1+1./R.^2).*sin(TH);

Ux=ur.*cos(TH)-uth.*sin(TH);
Uy=ur.*sin(TH)+uth.*cos(TH);

Umag=sqrt(Ux.^2+Uy.^2);

UmagS=Umag(:,1);
THS=TH(:,1);

[Umax,imax]=max(UmagS);
[Umin,imin]=min(UmagS);

thetaMax=THS(imax);
thetaMin=THS(imin);

figure
plot(THS,UmagS), grid on
hold on
plot(thetaMax,Umax,'ro',thetaMin,Umin,'bo')
title('Modulo de la velocidad sobre la frontera S')
xlabel('\theta'), ylabel('|\itu|')
hold off
}}

== Presión del fluido y ecuación de Bernoulli==

Se supone que la densidad del fluido es constante e igual a

<center><math>\rho=2</math></center>

y se cumple la ecuación de Bernoulli

<center><math>\frac12\rho|\vec u|^2+p=C</math></center>

Despejando la presión se obtiene

<center><math>p=C-\frac12\rho|\vec u|^2</math></center>

Tomando para la constante el valor <math>C=1</math>, la presión queda completamente determinada en toda la región del fluido.

El campo de presiones se ha representado en MATLAB mediante:

{{matlab|codigo=
rho=2;
C=1;

P=C-0.5*rho*Umag.^2;

PS=P(:,1);
[Pmax,ipmax]=max(PS);
[Pmin,ipmin]=min(PS);

thetaPmax=THS(ipmax);
thetaPmin=THS(ipmin);

figure
contourf(R.*cos(TH),R.*sin(TH),P,30,'LineStyle','none');
colorbar
hold on
viscircles([0 0],1,'Color','k')
axis equal
title('Presión p(x,y)')
xlabel('x'), ylabel('y')
hold off
}}
{| class="wikitable" style="border:0; margin:auto; background:none;"
| style="border:0;" | [[Archivo:Subir3.png|thumb|500px|Campo bidimensional del módulo de la velocidad ∣𝑢⃗(𝑥,𝑦)∣ alrededor del obstáculo circular]]
| style="border:0;" | [[Archivo:Subir2.png|thumb|500px|Módulo de la velocidad ∣𝑢⃗∣ sobre la frontera del obstáculo 𝑟=1 en función del ángulo 𝜃]]
|}
Se observa que la presión es máxima en los puntos de remanso, donde la velocidad es nula, y mínima en los puntos donde la velocidad alcanza su valor máximo. Esto concuerda con la ecuación de Bernoulli, ya que un aumento del módulo de la velocidad implica necesariamente una disminución de la presión.

==Velocidad y presión al rodear el obstáculo.==

Ya teniendo el campo de presiones se puede apreciar que tiene relación con la velocidad, ya que cuando disminuye la presión aumenta el módulo de la velocidad, tratándose de un resultado lógico debido a que cuando hablamos de un fluido incompresible su valor de Trinomio de Bernoulli ha de ser constante.
Lo podemos ver mediante la siguiente expresión: (Ec. Bernoulli).

<math>B=z+P/γ+v^2/2g</math>

B=valor del trinomio de Bernoulli

P= presión

z= altura

v=velocidad que lleva el fluido en determinado punto

g=valor de la gravedad

γ=valor del peso específico

Donde podemos ver que si la presión aumenta la velocidad tiene que disminuir y viceversa.

==Paradoja de D'Alembert.==

La magnitud pn representa la fuerza que el fluido ejerce de forma puntual sobre la superficie del obstáculo. Para estudiar si existe empuje neto en la dirección horizontal, basta analizar la suma de las proyecciones de todas esas fuerzas sobre el vector i. Si dicha suma resulta nula, concluiremos que no hay fuerza total en esa dirección.

Del apartado anterior conocemos ya la presión en la frontera del círculo unidad:
<math>p=10-∥u∥^2=9-4 sin⁡^2 θ-4sin⁡θ</math>

Cálculo del vector normal

Sobre la circunferencia de radio 1 el vector normal exterior es

<math>n=-eρ.</math>

Para expresarlo en coordenadas cartesianas utilizamos la matriz de cambio de base de polares a cartesianas. Tomando su traspuesta se obtiene:

\begin{pmatrix}
\cos \theta & -sen \theta & 0\\
sen \theta & \cos \theta & 0\\
0 & 0 & 1
\end{pmatrix}

Aplicando esta matriz a <math>(1,0,0)^T</math>, que corresponde al versor radial eρ, se obtiene:
\begin{pmatrix}
\cos \theta &\\
sen \theta &\\
0
\end{pmatrix}
<math>=cos⁡θi+sin⁡θj</math>

Al tomar el signo negativo para el normal exterior,

<math>n=cos⁡θi-sin⁡θj</math>,

de donde su proyección sobre i es:

<math>n⋅i=-cos⁡θ</math>.

Integración sobre la frontera

El borde del dominio se parametriza como
<math>(ρ,θ)=(1,θ), con θ∈[0,2π)</math>.

La contribución total de la presión en la dirección horizontal es:
<math>\int_{0}^{2π}p(n⋅i)dθ=\int_{0}^{2π}p(-cos⁡θ)dθ.</math>

Sustituyendo la expresión de p:
<math>\int_{0}^{2π}(9-4sin⁡^2 θ-4sin⁡θ)(-cos⁡θ)dθ.</math>

Integramos término a término:
<math>\int_{0}^{2π}-9cos⁡θdθ+\int_{0}^{2π}4sin⁡^2 θcos⁡θdθ+\int_{0}^{2π}4sin⁡θcos⁡θdθ.</math>

Cada una de estas integrales es periódica con una primitiva continua en todo el intervalo; evaluándolas entre 0 y 2π se obtiene:
<math>[-9sin⁡θ]_{0}^{2π}=0,[(4sin⁡^3 θ)/3]_{0}^{2π}=0,[2sin⁡^2 θ]_{0}^{2π}=0.</math>

Por tanto,
<math>\int_{0}^{2π}p(n⋅i)dθ=0.</math>

Conclusión

La componente horizontal de la fuerza total ejercida por el fluido sobre el obstáculo es nula. En consecuencia, el fluido no produce empuje neto en la dirección i.

==Cambio de función potencial==

En este apartado repetimos los procesos de los puntos (2), (3) y (4), pero tomando ahora como función potencial

<center><math>\varphi(\rho,\theta)=\left(\rho+\frac{1}{\rho}\right)\cos\theta+\frac{\theta}{4\pi}</math></center>

El obstáculo sigue siendo la circunferencia unidad centrada en el origen, y el fluido ocupa la región exterior, trabajándose de nuevo en coordenadas cilíndricas (polares).

===Representación de la función potencial y del campo de velocidades===

En primer lugar representamos la nueva función potencial <math>\varphi</math> sobre la región ocupada por el fluido, manteniendo el mismo mallado del apartado inicial. Esta representación permite visualizar cómo el término angular <math>\theta/(4\pi)</math> introduce una ligera inclinación en la superficie potencial respecto al caso anterior.

[[Archivo:Potencial9_g68.JPG|thumb|400px|right|Nueva función potencial en la región del fluido]]

La velocidad de las partículas del fluido viene dada de nuevo por el gradiente del potencial. En coordenadas cilíndricas, el campo de velocidades es

<center><math>\vec u(\rho,\theta,z)=\nabla\varphi
= u_\rho\,\vec e_\rho+u_\theta\,\vec e_\theta+u_z\,\vec e_z</math></center>

donde

<center><math>
u_\rho=\frac{\partial\varphi}{\partial\rho}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta, \qquad
u_\theta=\frac{1}{\rho}\frac{\partial\varphi}{\partial\theta}
= -\left(1+\frac{1}{\rho^2}\right)\sin\theta+\frac{1}{4\pi\rho},\qquad
u_z=0.
</math></center>

Para poder representarlo en MATLAB se pasa el campo a coordenadas cartesianas mediante la matriz de cambio de base de cilíndricas a cartesianas:

<center>
\begin{pmatrix}
\vec i\\
\vec j\\
\vec k
\end{pmatrix}
=
\begin{pmatrix}
\cos\theta & -\sin\theta & 0\\
\sin\theta & \cos\theta & 0\\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\vec e_\rho\\
\vec e_\theta\\
\vec e_z
\end{pmatrix}
</center>

Aplicando esta matriz al vector <math>(u_\rho,u_\theta,0)^T</math> se obtiene

<center><math>
\vec u(\vec i,\vec j,\vec k)
=\left[(1-\tfrac{1}{\rho^2})\cos^2\theta+(1+\tfrac{1}{\rho^2})\sin^2\theta
-\frac{\sin\theta}{4\pi\rho}\right]\vec i
+\left[(1-\tfrac{1}{\rho^2})\sin\theta\cos\theta
-(1+\tfrac{1}{\rho^2})\sin\theta\cos\theta
+\frac{\cos\theta}{4\pi\rho}\right]\vec j.
</math></center>

Sobre las curvas de nivel de <math>\varphi</math> se representa este nuevo campo de velocidades:

[[Archivo:CampoVelocidad9_g68.JPG|thumb|500px|right|Campo de velocidades sobre las líneas equipotenciales para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.1: Potencial y campo de velocidades para la nueva función φ

r = linspace(1,5,30); % Coordenada radial (rho)
a = linspace(0,2*pi,30); % Coordenada angular (theta)
[R,A] = meshgrid(r,a); % Mallado en (rho,theta)

% Conversión a coordenadas cartesianas
X = R.*cos(A);
Y = R.*sin(A);

% Nueva función potencial
phi = (R + 1./R).*cos(A) + A/(4*pi);

figure
surf(X,Y,phi);
shading interp
hold on
plot(1*cos(a),1*sin(a),'k','LineWidth',2); % Obstáculo
axis([-4 4 -4 4])
title('Nueva función potencial')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off

% Campo de velocidades (en cartesianas) sobre curvas de nivel de φ
figure
contour(X,Y,phi,30);
hold on

u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

quiver(X,Y,Cx,Cy)
plot(1*cos(a),1*sin(a),'k','LineWidth',1)
axis equal
axis([-4 4 -4 4])
title('Campo de velocidades (nuevo potencial)')
xlabel('EJE X'); ylabel('EJE Y');
colorbar
hold off
}}

Gracias a la representación anterior se observa que los vectores velocidad siguen siendo ortogonales a las líneas equipotenciales del nuevo potencial.

===Rotacional nulo y divergencia nula del nuevo campo===

Al ser <math>\vec u</math> el gradiente de un potencial escalar, su rotacional sigue siendo nulo:

<center><math>\nabla\times\vec u = \vec 0.</math></center>

En coordenadas cilíndricas se verifica explícitamente mediante la expresión

<center><math>
\nabla\times\vec u =
\frac{1}{\rho}
\begin{vmatrix}
\vec e_\rho & \vec e_\theta & \vec e_z\\
\frac{\partial}{\partial\rho} &
\frac{\partial}{\partial\theta} &
\frac{\partial}{\partial z}\\
u_\rho & \rho u_\theta & 0
\end{vmatrix}
= \vec 0,
</math></center>

lo que muestra que el flujo sigue siendo irrotacional.

Por otro lado, la divergencia del nuevo campo de velocidades viene dada por

<center><math>
\nabla\cdot\vec u(\rho,\theta,z)
=\frac{1}{\rho}\left\{
\frac{\partial}{\partial\rho}(\rho u_\rho)
+\frac{\partial}{\partial\theta}(u_\theta)
+\frac{\partial}{\partial z}(\rho u_z)
\right\},
</math></center>

y sustituyendo las expresiones de <math>u_\rho</math> y <math>u_\theta</math> se obtiene de nuevo

<center><math>\nabla\cdot\vec u = 0.</math></center>

La divergencia nula indica que el flujo sigue siendo incompresible: no hay creación ni destrucción de volumen de fluido en el dominio.

===Líneas de corriente del nuevo flujo===

Para determinar las líneas de corriente asociadas al nuevo campo de velocidades, consideramos el campo <math>\vec v</math> ortogonal a <math>\vec u</math>, definido por

<center><math>\vec v = \vec k \times \vec u.</math></center>

En coordenadas cilíndricas se tiene

<center>
<math>
\vec v
= \left(\left(1+\tfrac{1}{\rho^2}\right)\sin\theta - \tfrac{1}{4\pi\rho}\right)\vec e_\rho
+ \left(1-\tfrac{1}{\rho^2}\right)\cos\theta\,\vec e_\theta.
</math>
</center>

Dado que el rotacional de <math>\vec v</math> también es nulo, existe un potencial escalar <math>\psi</math> tal que

<center><math>\nabla\psi = \vec v.</math></center>

A partir de las componentes de <math>\vec v</math> se obtiene el siguiente sistema:

<center><math>
\frac{\partial\psi}{\partial\rho}
= \left(1+\frac{1}{\rho^2}\right)\sin\theta - \frac{1}{4\pi\rho},\qquad
\frac{1}{\rho}\frac{\partial\psi}{\partial\theta}
= \left(1-\frac{1}{\rho^2}\right)\cos\theta.
</math></center>

Integrando la primera ecuación respecto de <math>\rho</math> y comprobando la compatibilidad con la segunda, llegamos a

<center><math>
\psi(\rho,\theta)
= \left(\rho-\frac{1}{\rho}\right)\sin\theta
-\frac{1}{4\pi}\ln\rho + C,
</math></center>

donde <math>C</math> es una constante. Las curvas <math>\psi=\text{cte}</math> representan las líneas de corriente del nuevo flujo. Su representación numérica muestra cómo la introducción del término angular modifica el patrón de las corrientes, introduciendo una ligera asimetría respecto al caso inicial.

[[Archivo:LineasCorriente9_g68.JPG|thumb|400px|right|Líneas de corriente para el nuevo potencial]]

{{matlab|codigo=
% Apartado 9.3: Líneas de corriente del nuevo flujo

r = linspace(1,5,80);
a = linspace(0,2*pi,80);
[R,A] = meshgrid(r,a);

X = R.*cos(A);
Y = R.*sin(A);

% Campo de velocidades (reutilizamos las expresiones de 9.1)
u_rho = (1 - 1./R.^2).*cos(A);
u_theta = -(1 + 1./R.^2).*sin(A) + 1./(4*pi*R);

Cx = u_rho.*cos(A) - u_theta.*sin(A);
Cy = u_rho.*sin(A) + u_theta.*cos(A);

% Función de corriente psi
psi = (R - 1./R).*sin(A) - (1/(4*pi))*log(R);

figure
contourf(X,Y,psi,30,'LineStyle','-');
hold on
colorbar
quiver(X(1:3:end,1:3:end),Y(1:3:end,1:3:end), ...
Cx(1:3:end,1:3:end),Cy(1:3:end,1:3:end),1.3,'w');
plot(cos(a),sin(a),'k','LineWidth',2)
axis equal
axis([-4 4 -4 4])
title('Líneas de corriente para el nuevo potencial')
xlabel('EJE X'); ylabel('EJE Y');
hold off
}}

==Poster==
En el poster a continuación se muestra resumidamente el trabajo.
<center>[[Archivo:posterg68.JPG]]</center>

Archivo:Nuevafuncionpotencial.png

2025-12-07T22:52:13Z

Jeremy.garcia.herrera:

Archivo:Lineasdecorrienteparaelnuevopotencial.png

2025-12-07T22:52:01Z

Jeremy.garcia.herrera:

Archivo:Campodevelocidadesnuevopotencial.png

2025-12-07T22:51:37Z

Jeremy.garcia.herrera:

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:57:04Z

Jeremy.garcia.herrera: /* Velocidad extrema y puntos de remanso */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:32:42Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:32:10Z

Jeremy.garcia.herrera: /* Velocidad extrema y puntos de remanso */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:31:06Z

Jeremy.garcia.herrera: /* Velocidad extrema y puntos de remanso */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:30:09Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:27:55Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:26:15Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:08:43Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:07:12Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T21:05:26Z

Jeremy.garcia.herrera:

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:58:29Z

Jeremy.garcia.herrera: /* Velocidad extrema y puntos de remanso */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:57:59Z

Jeremy.garcia.herrera:

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:53:16Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:52:21Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:51:11Z

Jeremy.garcia.herrera: /* Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:50:22Z

Jeremy.garcia.herrera: /* 5. Velocidad extrema y puntos de remanso */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:35:17Z

Jeremy.garcia.herrera: /* 6. Presión del fluido y ecuación de Bernoulli */

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:27:31Z

Jeremy.garcia.herrera:

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:20:28Z

Jeremy.garcia.herrera:

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:17:46Z

Jeremy.garcia.herrera:

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-12-02T20:14:55Z

Jeremy.garcia.herrera:

Archivo:Subir1.png

2025-12-02T19:38:49Z

Jeremy.garcia.herrera:

Archivo:Subir2.png

2025-12-02T19:38:26Z

Jeremy.garcia.herrera:

Archivo:Subir3.png

2025-12-02T19:37:57Z

Jeremy.garcia.herrera:

Archivo:Imagen 3 jeremy.png

2025-11-30T09:06:52Z

Jeremy.garcia.herrera: Gráfica del módulo de la velocidad sobre la frontera del obstáculo

Gráfica del módulo de la velocidad sobre la frontera del obstáculo

Flujo alrededor de un obstáculo circular. (Grupo68)

2025-11-30T09:05:08Z

Jeremy.garcia.herrera:

Flujo de Poiseuille (Grupo 27-A)

2022-12-05T20:08:16Z

Jeremy.garcia.herrera: /* GRADIENTE DEL CAMPO DE TEMPERATURAS */

{{ TrabajoED | Flujo de Poiseuille (Grupo 27-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Pablo Díaz Paz-Albo Gonzalo de la Flor Fernández Jeremy García Herrera Juan Carlos Santos Expósito Emilio Villegas Maroto}}
{{ Beta }}
==INTRODUCCIÓN==
La ley de Poiseuille o ley de Hagen-Poiseuille es una ley que permite determinar el flujo laminar estacionario de un líquido incompresible y uniformemente viscoso (fluido newtoniano) a través de un tubo cilíndrico de sección circular constante.

Para la realización de este trabajo hemos utilizado el programa informático Matlab, para poder representar los resultados de una manera visual y ayudar al usuario a comprender mejor las interpretaciones.

==CAMPO DE VELOCIDADES==
Para la representación del campo se han tomado los siguientes valores:

<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>

Sustituyendo en la expresión del campo obtenemos:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>
Para su representación se ha introducido en Matlab, el siguiente código:
{{matlab|codigo=
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
ux=(16*xx.^2)/3-6.*xx;
uy=0.*yy;
hold on
quiver(xx,yy,ux,uy)
axis([0,5,0,12])
hold off
view(2)
title('CAMPO DE VELOCIDADES')
}}
[[Archivo:Campo_de_velocidades_poiseuille.png|400px|miniaturadeimagen|centro|Campo de velocidades]]

=== CÁLCULO DE LA VELOCIDAD MÁXIMA DEL FLUIDO===
Los puntos de velocidad máxima se obtienen igualando a 0, la primera derivada parcial del campo:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>

donde:

<center><math> \frac{\partial \vec{u}} {\partial \rho}=2\rho(\rho-2)\vec{e_z} </math></center>
<center><math> \frac{\partial \vec{u} }{\partial \theta}=0 </math></center>
<center><math> \frac{\partial \vec{u} }{\partial z}=0 </math></center>

Si igualamos a 0, <math> 2\rho(\rho-2)=0</math> obtenemos que los puntos donde se produce la máxima velocidad son <math> \rho =0 \ {y} \ \rho = 2 </math>

==ROTACIONAL==
Para la realización del cálculo del rotacional del campo, utilizaremos la siguiente fórmula:

<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}</math></big></center>

Sustituyendo en la fórmula:
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & -(\frac{2\rho^3}{3}-\rho^2)\frac{\cdot(p_1-p_2)}{μ} \end{vmatrix}</math></big></center>
Obtenemos que <math>\nabla\times\vec{u}=(\frac{8\rho^2}{3}-3\rho)\frac{(p_1-p_2)}{μ}\vec{e_\theta}</math>, y dando los siguientes valores:
<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>
Tenemos que <math>\nabla\times\vec u=(\frac{16\rho^2}{3}-6\rho)\vec{e_\theta}</math>, siendo el módulo del rotacional, <math>\left | \nabla\times\vec u \right |= \frac{16\rho^2}{3}-6\rho</math>, como se observa en la expresión, depende del valor de <math>\rho</math>. Tomando los valores mínimos en los extremos del intervalo en el que está definido <math>\rho</math>, <math> [0,1]</math>.

Para visualizar esto de manera gráfica, nuevamente a Octave,
{{matlab|codigo=
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs((16*xx.^2)/3-6.*xx);
surf(xx,yy,rot)
colorbar
view(2)
axis([0,5,0,12])
}}

[[Archivo:Rotacionaldepouseuille.png|400px|miniaturadeimagen|centro|Rotacional del campo velocidad]]

== CAMPO DE TEMPERATURAS==

La temperatura a la que está sometido el fluido viene determinada por el campo escalar:
<center><math> T(ρ{,}θ{,}z)=1+(ρ-1/2)^2 e^{-(z-1)^2}</math></center>

Como se observa, la expresión del campo está en coordenadas porlares. Teniendo en cuenta que las relación para pasar de coordenadas polares a cartesianas es:
: <math>x=ρcosθ</math>
: <math>y=ρsenθ</math>

Pero como nosotros estamos en los ejes <math>\{\vec{j} {,} \vec{k} \}</math>, usaremos:

: <math>y=ρcosθ</math>
: <math>z=ρsenθ</math>
Podemos expresar el campo en coordenadas cartesianas, que es el sistema en el cual estamos trabajando, resultando:
<center><math> T(y{,}z)=1+(y^2+z^2-1/2)^2 e^{-(z-1)^2}</math></center>

''' REPRESENTACIÓN DEL CAMPO DE TEMPERATURAS '''

El campo de temperaturas es un campo escalar que representa la distribución espacial de la temperatura asociando un valor a cada punto del espacio. Depende de la posición del punto, y del tiempo (t). En nuestro caso, el campo de temperaturas como ya se a expuesto viene dado por la expresión <math> T(y{,}z)=1+(y^2+z^2-1/2) e^{-((z-1)^2)}</math>, por lo tanto no depende del tiempo, se dice que es estacionario, sino exclusivamente de las componentes espaciales <math>(y{,}z)</math>

Para la representación se ha implementado el siguiente código en Octave:
{{matlab|codigo=
y=0:0.05:8;
z=-2:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
pcolor(Y,Z,p);
shading flat
grid on
axis([0,8,-1,10]);
colorbar
figure(2)
contour(Y,Z,p,10,'k');
grid on
axis([0,8,-1,10]);
}}
[[Archivo:Temperatura27-1.PNG|400px|miniaturadeimagen|centro| Campo de temperaturas del fluido]]
<center>'''INTERPRETACIÓN''' </center>
En el gráfico se puede observar que la temperatura es mayor cuándo se representa con colores cálidos, esto sucede en los valores próximos a <math>z=0</math> y <math>z=1</math>.

''' REPRESENTACIÓN DE LAS CURVAS DE NIVEL '''
[[Archivo:Temperatura27-2.PNG|400px|miniaturadeimagen|centro| Curvas de nivel de campo]]

El código utilizado para la representación de este gráfico se incluyó con el de la representación del campo de temperaturas.

<center>'''INTERPRETACIÓN''' </center>

En el gráfico aparecen representadas las curvas de nivel, las curvas de nivel varían de manera geométrica estando más próximas entre sí cuando la temperatura aumenta, en torono a los puntos <math>z=1</math> y más alejadas cuando la temperatura disminuya.

La representaciónes anteriores permite determinar que en el punto <math>(y,z)= (8,1)</math> la temperatura es máxima.

===GRADIENTE DEL CAMPO DE TEMPERATURAS===
El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente <math>\nabla T</math>, en nuestro caso:

<center><math>\nabla T(y,z)=(4y^3e^{-(z-1)^2}+4yz^2e^{-(z-1)^2}-2ye^{-(z-1)^2})\vec{j}</math> <math>+(4z^3e^{-(z-1)^2}-2z^4e^{-(z-1)^2}+4zy^2e^{-(z-1)^2}-4y^2z^2e^{-(z-1)^2} + 2ze^{-(z-1)^2}-2ze^{-(z-1)^2}) \vec{k}</math></center>

Como observación, el gradiente analíticamente se calcula como <math>\frac{\partial T(y,z)}{\partial y}</math> siendo ésta la compontente <math>\vec{j}</math> del campo vectorial gradiente y <math>\frac{\partial T(y,z)}{\partial z }</math> la componente <math>\vec{k}</math>.

''' REPRESENTACIÓN GRÁFICA DEL GRADIENTE''''''

El gráfico del gradiente de la temperatura, se ha representado utilizando el siguiente código:
{{Matlab|codigo=
y=0:0.05:8;
z=0:0.05:1;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+((Z.^2).*exp(-((Y.^2)+(Z.^2)-1/2).^2));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
axis([0,8,-1,2]);
shading flat
grid on
hold off
}}
[[Archivo:Temperatura27-4.PNG|400px|miniaturadeimagen|centro| Campo vectorial gradiente]]
Código Octave utilizado:
{{Matlab|codigo=

y=0:0.05:8;
z=0:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
contour(Y,Z,p,'k');
axis([0,8,-1,10]);
shading flat
grid on
hold off

}}
[[Archivo:Temperatura27-5.PNG|400px|thumb|centro|Gradiente ortogonal a las curvas de nivel]]

Con este gráfico afirmamos que las curvas de nivel son ortogonales al gradiente de la temperatura, es el mismo gráfico que el anterior pero más ampliado para poder observar mejor la ortogonalidad.

<center>'''INTERPRETACIÓN''' </center>

Para estudiar la variación de temperatura a lo largo del canal analizamos el gradiente que indica la direccion de crecimiento de la función temperatura en cada uno de los puntos. Como ya se ha mencionado antes, junto al gradiente aparecen las curvas de nivel de la temperatura. Se puede comprobar que en torno al punto <math>z=1</math> el módulo del gradiente aumenta, lo que indica que la temperatura tendrá una variación más rápida en esta zona. Se aprecia claramente como el gradiente es ortogonal a las curvas de nivel, esto se debe a que el gradiente indica la dirección de crecimiento y sentido de los valores del campo en cada punto y las curvas de nivel son el lugar geométrico de los puntos equipotenciales.

==CAUDAL==
El caudal es el volumen de un fluido que pasa a través de la tubería estudiada por unidad de tiempo y se calcula mediante la integral de superficie siguiente:

<center><math>\int_{S}\vec{v} \cdot d\vec{S}</math></center>

donde <math>\vec{v}</math> es el campo de velocidades. Al tratarse de una integral de superficie hay que tener en cuenta el vector normal que es perpendicular a la superficie.

En nuestro caso, el campo de velocidades es:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>,

Esa expresión resulta de sustituir los valores siguientes en la expresión incial del campo:

<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>

La profundidad del caudal es de 1 metro, por lo que al parametrizar nos da:
: <math>\rho= 15 </math>
: <math>\theta=u</math><math>[0,2\pi]</math>
: <math>z=v </math> <math>[0,\inf]</math>
Por lo que:
: <math>ru=(0,1,0)</math>
: <math>rv=(0,0,1)</math>
: <math>|ru\times rv|=(1,0,0)</math>

'''CÁLCULO DEL CAUDAL'''

<center><math>\int_{S}\vec{v} \cdot d\vec{S}=\int_{S}\vec{u} \cdot |ru\times rv|=\int_{0}^{\inf} \int_{0}^{2\pi}\vec{u}(\rho,\theta,z)\cdot |\vec{e_\theta}\times \vec{e_z}|du dv=\int_{0}^{\inf} \int_{0}^{2\pi}-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}(\vec{e_\rho})=0\frac{m^3}{s}</math></center>

[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]

Flujo de Poiseuille (Grupo 27-A)

2022-12-05T19:56:30Z

Jeremy.garcia.herrera: /* CAMPO DE TEMPERATURAS */

{{ TrabajoED | Flujo de Poiseuille (Grupo 27-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Pablo Díaz Paz-Albo Gonzalo de la Flor Fernández Jeremy García Herrera Juan Carlos Santos Expósito Emilio Villegas Maroto}}
{{ Beta }}
==INTRODUCCIÓN==
La ley de Poiseuille o ley de Hagen-Poiseuille es una ley que permite determinar el flujo laminar estacionario de un líquido incompresible y uniformemente viscoso (fluido newtoniano) a través de un tubo cilíndrico de sección circular constante.

Para la realización de este trabajo hemos utilizado el programa informático Matlab, para poder representar los resultados de una manera visual y ayudar al usuario a comprender mejor las interpretaciones.

==CAMPO DE VELOCIDADES==
Para la representación del campo se han tomado los siguientes valores:

<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>

Sustituyendo en la expresión del campo obtenemos:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>
Para su representación se ha introducido en Matlab, el siguiente código:
{{matlab|codigo=
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
ux=(16*xx.^2)/3-6.*xx;
uy=0.*yy;
hold on
quiver(xx,yy,ux,uy)
axis([0,5,0,12])
hold off
view(2)
title('CAMPO DE VELOCIDADES')
}}
[[Archivo:Campo_de_velocidades_poiseuille.png|400px|miniaturadeimagen|centro|Campo de velocidades]]

=== CÁLCULO DE LA VELOCIDAD MÁXIMA DEL FLUIDO===
Los puntos de velocidad máxima se obtienen igualando a 0, la primera derivada parcial del campo:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>

donde:

<center><math> \frac{\partial \vec{u}} {\partial \rho}=2\rho(\rho-2)\vec{e_z} </math></center>
<center><math> \frac{\partial \vec{u} }{\partial \theta}=0 </math></center>
<center><math> \frac{\partial \vec{u} }{\partial z}=0 </math></center>

Si igualamos a 0, <math> 2\rho(\rho-2)=0</math> obtenemos que los puntos donde se produce la máxima velocidad son <math> \rho =0 \ {y} \ \rho = 2 </math>

==ROTACIONAL==
Para la realización del cálculo del rotacional del campo, utilizaremos la siguiente fórmula:

<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}</math></big></center>

Sustituyendo en la fórmula:
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & -(\frac{2\rho^3}{3}-\rho^2)\frac{\cdot(p_1-p_2)}{μ} \end{vmatrix}</math></big></center>
Obtenemos que <math>\nabla\times\vec{u}=(\frac{8\rho^2}{3}-3\rho)\frac{(p_1-p_2)}{μ}\vec{e_\theta}</math>, y dando los siguientes valores:
<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>
Tenemos que <math>\nabla\times\vec u=(\frac{16\rho^2}{3}-6\rho)\vec{e_\theta}</math>, siendo el módulo del rotacional, <math>\left | \nabla\times\vec u \right |= \frac{16\rho^2}{3}-6\rho</math>, como se observa en la expresión, depende del valor de <math>\rho</math>. Tomando los valores mínimos en los extremos del intervalo en el que está definido <math>\rho</math>, <math> [0,1]</math>.

Para visualizar esto de manera gráfica, nuevamente a Octave,
{{matlab|codigo=
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs((16*xx.^2)/3-6.*xx);
surf(xx,yy,rot)
colorbar
view(2)
axis([0,5,0,12])
}}

[[Archivo:Rotacionaldepouseuille.png|400px|miniaturadeimagen|centro|Rotacional del campo velocidad]]

== CAMPO DE TEMPERATURAS==

La temperatura a la que está sometido el fluido viene determinada por el campo escalar:
<center><math> T(ρ{,}θ{,}z)=1+(ρ-1/2)^2 e^{-(z-1)^2}</math></center>

Como se observa, la expresión del campo está en coordenadas porlares. Teniendo en cuenta que las relación para pasar de coordenadas polares a cartesianas es:
: <math>x=ρcosθ</math>
: <math>y=ρsenθ</math>

Pero como nosotros estamos en los ejes <math>\{\vec{j} {,} \vec{k} \}</math>, usaremos:

: <math>y=ρcosθ</math>
: <math>z=ρsenθ</math>
Podemos expresar el campo en coordenadas cartesianas, que es el sistema en el cual estamos trabajando, resultando:
<center><math> T(y{,}z)=1+(y^2+z^2-1/2)^2 e^{-(z-1)^2}</math></center>

''' REPRESENTACIÓN DEL CAMPO DE TEMPERATURAS '''

El campo de temperaturas es un campo escalar que representa la distribución espacial de la temperatura asociando un valor a cada punto del espacio. Depende de la posición del punto, y del tiempo (t). En nuestro caso, el campo de temperaturas como ya se a expuesto viene dado por la expresión <math> T(y{,}z)=1+(y^2+z^2-1/2) e^{-((z-1)^2)}</math>, por lo tanto no depende del tiempo, se dice que es estacionario, sino exclusivamente de las componentes espaciales <math>(y{,}z)</math>

Para la representación se ha implementado el siguiente código en Octave:
{{matlab|codigo=
y=0:0.05:8;
z=-2:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
pcolor(Y,Z,p);
shading flat
grid on
axis([0,8,-1,10]);
colorbar
figure(2)
contour(Y,Z,p,10,'k');
grid on
axis([0,8,-1,10]);
}}
[[Archivo:Temperatura27-1.PNG|400px|miniaturadeimagen|centro| Campo de temperaturas del fluido]]
<center>'''INTERPRETACIÓN''' </center>
En el gráfico se puede observar que la temperatura es mayor cuándo se representa con colores cálidos, esto sucede en los valores próximos a <math>z=0</math> y <math>z=1</math>.

''' REPRESENTACIÓN DE LAS CURVAS DE NIVEL '''
[[Archivo:Temperatura27-2.PNG|400px|miniaturadeimagen|centro| Curvas de nivel de campo]]

El código utilizado para la representación de este gráfico se incluyó con el de la representación del campo de temperaturas.

<center>'''INTERPRETACIÓN''' </center>

En el gráfico aparecen representadas las curvas de nivel, las curvas de nivel varían de manera geométrica estando más próximas entre sí cuando la temperatura aumenta, en torono a los puntos <math>z=1</math> y más alejadas cuando la temperatura disminuya.

La representaciónes anteriores permite determinar que en el punto <math>(y,z)= (8,1)</math> la temperatura es máxima.

===GRADIENTE DEL CAMPO DE TEMPERATURAS===
El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente <math>\nabla T</math>, en nuestro caso:

<center><math>\nabla T(y,z)=(4y^3e^[-(z-1)^2]+4yz^2e^(-(z-1)^2)-2ye^(-(z-1)^2))vec{j} + (4z^3e^(-(z-1)^2)-2z^4e^(-(z-1)^2)+4zy^2e^(-(z-1)^2)-4y^2z^2e^(-(z-1)^2) + 2ze^(-(z-1)^2)-2ze^(-(z-1)^2))vec{k} </math></center>

Como observación, el gradiente analíticamente se calcula como <math>\frac{\partial T(y,z)}{\partial y}</math> siendo ésta la compontente <math>\vec{j}</math> del campo vectorial gradiente y <math>\frac{\partial T(y,z)}{\partial z }</math> la componente <math>\vec{k}</math>.

''' REPRESENTACIÓN GRÁFICA DEL GRADIENTE''''''

El gráfico del gradiente de la temperatura, se ha representado utilizando el siguiente código:
{{Matlab|codigo=
y=0:0.05:8;
z=0:0.05:1;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+((Z.^2).*exp(-((Y.^2)+(Z.^2)-1/2).^2));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
axis([0,8,-1,2]);
shading flat
grid on
hold off
}}
[[Archivo:Temperatura27-4.PNG|400px|miniaturadeimagen|centro| Campo vectorial gradiente]]
Código Octave utilizado:
{{Matlab|codigo=

y=0:0.05:8;
z=0:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
contour(Y,Z,p,'k');
axis([0,8,-1,10]);
shading flat
grid on
hold off

}}
[[Archivo:Temperatura27-5.PNG|400px|thumb|centro|Gradiente ortogonal a las curvas de nivel]]

Con este gráfico afirmamos que las curvas de nivel son ortogonales al gradiente de la temperatura, es el mismo gráfico que el anterior pero más ampliado para poder observar mejor la ortogonalidad.

<center>'''INTERPRETACIÓN''' </center>

Para estudiar la variación de temperatura a lo largo del canal analizamos el gradiente que indica la direccion de crecimiento de la función temperatura en cada uno de los puntos. Como ya se ha mencionado antes, junto al gradiente aparecen las curvas de nivel de la temperatura. Se puede comprobar que en torno al punto <math>z=1</math> el módulo del gradiente aumenta, lo que indica que la temperatura tendrá una variación más rápida en esta zona. Se aprecia claramente como el gradiente es ortogonal a las curvas de nivel, esto se debe a que el gradiente indica la dirección de crecimiento y sentido de los valores del campo en cada punto y las curvas de nivel son el lugar geométrico de los puntos equipotenciales.

==CAUDAL==
El caudal es el volumen de un fluido que pasa a través de la tubería estudiada por unidad de tiempo y se calcula mediante la integral de superficie siguiente:

<center><math>\int_{S}\vec{v} \cdot d\vec{S}</math></center>

donde <math>\vec{v}</math> es el campo de velocidades. Al tratarse de una integral de superficie hay que tener en cuenta el vector normal que es perpendicular a la superficie.

En nuestro caso, el campo de velocidades es:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>,

Esa expresión resulta de sustituir los valores siguientes en la expresión incial del campo:

<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>

La profundidad del caudal es de 1 metro, por lo que al parametrizar nos da:
: <math>\rho= 15 </math>
: <math>\theta=u</math><math>[0,2\pi]</math>
: <math>z=v </math> <math>[0,\inf]</math>
Por lo que:
: <math>ru=(0,1,0)</math>
: <math>rv=(0,0,1)</math>
: <math>|ru\times rv|=(1,0,0)</math>

'''CÁLCULO DEL CAUDAL'''

<center><math>\int_{S}\vec{v} \cdot d\vec{S}=\int_{S}\vec{u} \cdot |ru\times rv|=\int_{0}^{\inf} \int_{0}^{2\pi}\vec{u}(\rho,\theta,z)\cdot |\vec{e_\theta}\times \vec{e_z}|du dv=\int_{0}^{\inf} \int_{0}^{2\pi}-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}(\vec{e_\rho})=0\frac{m^3}{s}</math></center>

[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]

Flujo de Poiseuille (Grupo 27-A)

2022-12-05T19:53:43Z

Jeremy.garcia.herrera: /* GRADIENTE DEL CAMPO DE TEMPERATURAS */

{{ TrabajoED | Flujo de Poiseuille (Grupo 27-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Pablo Díaz Paz-Albo Gonzalo de la Flor Fernández Jeremy García Herrera Juan Carlos Santos Expósito Emilio Villegas Maroto}}
{{ Beta }}
==INTRODUCCIÓN==
La ley de Poiseuille o ley de Hagen-Poiseuille es una ley que permite determinar el flujo laminar estacionario de un líquido incompresible y uniformemente viscoso (fluido newtoniano) a través de un tubo cilíndrico de sección circular constante.

Para la realización de este trabajo hemos utilizado el programa informático Matlab, para poder representar los resultados de una manera visual y ayudar al usuario a comprender mejor las interpretaciones.

==CAMPO DE VELOCIDADES==
Para la representación del campo se han tomado los siguientes valores:

<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>

Sustituyendo en la expresión del campo obtenemos:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>
Para su representación se ha introducido en Matlab, el siguiente código:
{{matlab|codigo=
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
ux=(16*xx.^2)/3-6.*xx;
uy=0.*yy;
hold on
quiver(xx,yy,ux,uy)
axis([0,5,0,12])
hold off
view(2)
title('CAMPO DE VELOCIDADES')
}}
[[Archivo:Campo_de_velocidades_poiseuille.png|400px|miniaturadeimagen|centro|Campo de velocidades]]

=== CÁLCULO DE LA VELOCIDAD MÁXIMA DEL FLUIDO===
Los puntos de velocidad máxima se obtienen igualando a 0, la primera derivada parcial del campo:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>

donde:

<center><math> \frac{\partial \vec{u}} {\partial \rho}=2\rho(\rho-2)\vec{e_z} </math></center>
<center><math> \frac{\partial \vec{u} }{\partial \theta}=0 </math></center>
<center><math> \frac{\partial \vec{u} }{\partial z}=0 </math></center>

Si igualamos a 0, <math> 2\rho(\rho-2)=0</math> obtenemos que los puntos donde se produce la máxima velocidad son <math> \rho =0 \ {y} \ \rho = 2 </math>

==ROTACIONAL==
Para la realización del cálculo del rotacional del campo, utilizaremos la siguiente fórmula:

<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ u_\rho & \rho u_\theta & u_z \end{vmatrix}</math></big></center>

Sustituyendo en la fórmula:
<center> <big> <math>\nabla\times\vec u= \frac{1}{\rho}\begin{vmatrix} \vec{e_\rho} & \rho \vec{e_\theta} & \vec{e_z} \\ \frac{ \partial}{\partial \rho} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z}\\ 0 & 0 & -(\frac{2\rho^3}{3}-\rho^2)\frac{\cdot(p_1-p_2)}{μ} \end{vmatrix}</math></big></center>
Obtenemos que <math>\nabla\times\vec{u}=(\frac{8\rho^2}{3}-3\rho)\frac{(p_1-p_2)}{μ}\vec{e_\theta}</math>, y dando los siguientes valores:
<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>
Tenemos que <math>\nabla\times\vec u=(\frac{16\rho^2}{3}-6\rho)\vec{e_\theta}</math>, siendo el módulo del rotacional, <math>\left | \nabla\times\vec u \right |= \frac{16\rho^2}{3}-6\rho</math>, como se observa en la expresión, depende del valor de <math>\rho</math>. Tomando los valores mínimos en los extremos del intervalo en el que está definido <math>\rho</math>, <math> [0,1]</math>.

Para visualizar esto de manera gráfica, nuevamente a Octave,
{{matlab|codigo=
x=0:0.1:2;
y=-1:0.1:10;
[xx,yy]=meshgrid(x,y);
rot=abs((16*xx.^2)/3-6.*xx);
surf(xx,yy,rot)
colorbar
view(2)
axis([0,5,0,12])
}}

[[Archivo:Rotacionaldepouseuille.png|400px|miniaturadeimagen|centro|Rotacional del campo velocidad]]

== CAMPO DE TEMPERATURAS==

La temperatura a la que está sometido el fluido viene determinada por el campo escalar:
<center><math> T(ρ{,}θ{,}z)=1+(ρ-1/2)^2 e^{-(z-1)^2}</math></center>

Como se observa, la expresión del campo está en coordenadas porlares. Teniendo en cuenta que las relación para pasar de coordenadas polares a cartesianas es:
: <math>x=ρcosθ</math>
: <math>y=ρsenθ</math>

Pero como nosotros estamos en los ejes <math>\{\vec{j} {,} \vec{k} \}</math>, usaremos:

: <math>y=ρcosθ</math>
: <math>z=ρsenθ</math>
Podemos expresar el campo en coordenadas cartesianas, que es el sistema en el cual estamos trabajando, resultando:
<center><math> T(y{,}z)=1+(y^2+z^2-1/2) e^{-(z-1)}</math></center>

''' REPRESENTACIÓN DEL CAMPO DE TEMPERATURAS '''

El campo de temperaturas es un campo escalar que representa la distribución espacial de la temperatura asociando un valor a cada punto del espacio. Depende de la posición del punto, y del tiempo (t). En nuestro caso, el campo de temperaturas como ya se a expuesto viene dado por la expresión <math> T(y{,}z)=1+(y^2+z^2-1/2) e^{-((z-1)^2)}</math>, por lo tanto no depende del tiempo, se dice que es estacionario, sino exclusivamente de las componentes espaciales <math>(y{,}z)</math>

Para la representación se ha implementado el siguiente código en Octave:
{{matlab|codigo=
y=0:0.05:8;
z=-2:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
pcolor(Y,Z,p);
shading flat
grid on
axis([0,8,-1,10]);
colorbar
figure(2)
contour(Y,Z,p,10,'k');
grid on
axis([0,8,-1,10]);
}}
[[Archivo:Temperatura27-1.PNG|400px|miniaturadeimagen|centro| Campo de temperaturas del fluido]]
<center>'''INTERPRETACIÓN''' </center>
En el gráfico se puede observar que la temperatura es mayor cuándo se representa con colores cálidos, esto sucede en los valores próximos a <math>z=0</math> y <math>z=1</math>.

''' REPRESENTACIÓN DE LAS CURVAS DE NIVEL '''
[[Archivo:Temperatura27-2.PNG|400px|miniaturadeimagen|centro| Curvas de nivel de campo]]

El código utilizado para la representación de este gráfico se incluyó con el de la representación del campo de temperaturas.

<center>'''INTERPRETACIÓN''' </center>

En el gráfico aparecen representadas las curvas de nivel, las curvas de nivel varían de manera geométrica estando más próximas entre sí cuando la temperatura aumenta, en torono a los puntos <math>z=1</math> y más alejadas cuando la temperatura disminuya.

La representaciónes anteriores permite determinar que en el punto <math>(y,z)= (8,1)</math> la temperatura es máxima.

===GRADIENTE DEL CAMPO DE TEMPERATURAS===
El análisis de la variación de la temperatura a lo largo del canal se realiza analizando el gradiente <math>\nabla T</math>, en nuestro caso:

<center><math>\nabla T(y,z)=(4y^3e^[-(z-1)^2]+4yz^2e^(-(z-1)^2)-2ye^(-(z-1)^2))vec{j} + (4z^3e^(-(z-1)^2)-2z^4e^(-(z-1)^2)+4zy^2e^(-(z-1)^2)-4y^2z^2e^(-(z-1)^2) + 2ze^(-(z-1)^2)-2ze^(-(z-1)^2))vec{k} </math></center>

Como observación, el gradiente analíticamente se calcula como <math>\frac{\partial T(y,z)}{\partial y}</math> siendo ésta la compontente <math>\vec{j}</math> del campo vectorial gradiente y <math>\frac{\partial T(y,z)}{\partial z }</math> la componente <math>\vec{k}</math>.

''' REPRESENTACIÓN GRÁFICA DEL GRADIENTE''''''

El gráfico del gradiente de la temperatura, se ha representado utilizando el siguiente código:
{{Matlab|codigo=
y=0:0.05:8;
z=0:0.05:1;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+((Z.^2).*exp(-((Y.^2)+(Z.^2)-1/2).^2));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
axis([0,8,-1,2]);
shading flat
grid on
hold off
}}
[[Archivo:Temperatura27-4.PNG|400px|miniaturadeimagen|centro| Campo vectorial gradiente]]
Código Octave utilizado:
{{Matlab|codigo=

y=0:0.05:8;
z=0:0.05:10;
[Y,Z]=meshgrid(y,z);
figure (1)
p=1+(((Y.^2+Z.^2-(1/2)).^2).*exp(-((Z-1).^2)));
[pY,pZ]=gradient(p);
hold on
quiver(Y,Z,pY,pZ)
contour(Y,Z,p,'k');
axis([0,8,-1,10]);
shading flat
grid on
hold off

}}
[[Archivo:Temperatura27-5.PNG|400px|thumb|centro|Gradiente ortogonal a las curvas de nivel]]

Con este gráfico afirmamos que las curvas de nivel son ortogonales al gradiente de la temperatura, es el mismo gráfico que el anterior pero más ampliado para poder observar mejor la ortogonalidad.

<center>'''INTERPRETACIÓN''' </center>

Para estudiar la variación de temperatura a lo largo del canal analizamos el gradiente que indica la direccion de crecimiento de la función temperatura en cada uno de los puntos. Como ya se ha mencionado antes, junto al gradiente aparecen las curvas de nivel de la temperatura. Se puede comprobar que en torno al punto <math>z=1</math> el módulo del gradiente aumenta, lo que indica que la temperatura tendrá una variación más rápida en esta zona. Se aprecia claramente como el gradiente es ortogonal a las curvas de nivel, esto se debe a que el gradiente indica la dirección de crecimiento y sentido de los valores del campo en cada punto y las curvas de nivel son el lugar geométrico de los puntos equipotenciales.

==CAUDAL==
El caudal es el volumen de un fluido que pasa a través de la tubería estudiada por unidad de tiempo y se calcula mediante la integral de superficie siguiente:

<center><math>\int_{S}\vec{v} \cdot d\vec{S}</math></center>

donde <math>\vec{v}</math> es el campo de velocidades. Al tratarse de una integral de superficie hay que tener en cuenta el vector normal que es perpendicular a la superficie.

En nuestro caso, el campo de velocidades es:

<center><math>\vec{u}(\rho,\theta,z)=-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}</math></center>,

Esa expresión resulta de sustituir los valores siguientes en la expresión incial del campo:

<center><math>p_1=2 \ {,} \ p_2=1 \ {,} \ μ=1</math></center>

La profundidad del caudal es de 1 metro, por lo que al parametrizar nos da:
: <math>\rho= 15 </math>
: <math>\theta=u</math><math>[0,2\pi]</math>
: <math>z=v </math> <math>[0,\inf]</math>
Por lo que:
: <math>ru=(0,1,0)</math>
: <math>rv=(0,0,1)</math>
: <math>|ru\times rv|=(1,0,0)</math>

'''CÁLCULO DEL CAUDAL'''

<center><math>\int_{S}\vec{v} \cdot d\vec{S}=\int_{S}\vec{u} \cdot |ru\times rv|=\int_{0}^{\inf} \int_{0}^{2\pi}\vec{u}(\rho,\theta,z)\cdot |\vec{e_\theta}\times \vec{e_z}|du dv=\int_{0}^{\inf} \int_{0}^{2\pi}-(\frac{2\rho^3}{3}-\rho^2)\vec{e_z}(\vec{e_\rho})=0\frac{m^3}{s}</math></center>

[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]