Flujo de Couette entre dos tubos concéntricos (Grupo 25)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Flujo de Couette entre dos tubos concéntricos . Grupo 25 |
| Asignatura | Teoría de Campos |
| Curso | 2025-26 |
| Autores | Andrés Gaitán Tsukasov, Pablo Casado del Campo, Esteban Gabaldón Hermoso , Carlos Ascanio Martín, Aarón Pérez Luna |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Flujo de Couette entre dos tubos concéntricos
Introducción
El flujo de Couette entre dos tubos concentricos, designa el movimiento de un fluido viscoso entre dos cilindros coaxiales de diferente radio, siendo dominado así el movimiento por la rotación de los cilindros alrederor de una eje común. En ingenieria, el flujo de Couette es sumamente importante ya que es uno de los pocos fluidos viscosos de los que se conoce la solución exacta de la ecuacion de Navier-Stokes, volviendolo fundamental en el estudio de propiedades de los fluidos y estabilidad de los flujos.
En este trabajo partiremos de un flujo generado por la rotación en sentido antihorario del cilindro interior con velocidad angular constante [math]\vec{w}_i[/math] y la rotación en sentido horario del cilindro exterior con velocidad angular [math]\vec{w}_e[/math]. Ambos cilindros tienen su eje [math]OX_3[/math] y sus radios [math]\rho[/math]=2 para el cilindro exterior y [math]\rho[/math]=1 para el interior. A partir de ello, con el campo de velocidades entre ambos cilindros, estudiaremos el flujo de temperaturas resultante, el gradiente y el caudal, entre otras.
Contenido
- 1 Mallado de la sección transversal
- 2 Calculo de velocidades. Ecuación de Navier-Stokes
- 2.1 Velocidad de las partículas del fluido y particularización de la ecuación
- 2.2 Cálculo del gradiente de una divergencia
- 2.3 Cálculo del rotacional del rotacional
- 2.4 Solución de la ecuacion de Navier-Stokes de forma diferencial
- 2.5 Comprobación de la ecuación obtenida
- 2.6 Calculo de los parámetros a y b
- 2.7 Condición de incompresibilidad
- 3 Representación del campo de velocidades
- 4 Líneas de Corriente
- 5 Gradiente de la temperatura
- 6 Caudal de una sección longitudinal
1 Mallado de la sección transversal
En primer lugar, se toma la sección de ambos cilindros en el plano [math]x_3[/math]=0 y se representa en matlab de acuerdo al parámetro angular [math]\theta \in [0,2\pi][/math] y a los parámetros cartesianos [math](x,y) \in [-2,2] \times [-2,2][/math].
% MALLADO DE LA SECCIÓN TRANSVERSAL
figure;
rho=1:0.1:2; %SE ESTABLECEN PARÁMETROS
theta_malla = linspace(0, 2*pi, 50);
[RHO,THETA]=meshgrid(rho,theta_malla); %OBTENCION DE LA MALLA
rho_int=1;
rho_ext=2;
hold on
grid on
x=RHO.*cos(THETA);
y=RHO.*sin(THETA);
mesh(x,y,0*y,'EdgeColor','g');
x_ext = rho_ext.*cos(theta_malla); %OBTENCIÓN DE COORDENADAS
y_ext = rho_ext.*sin(theta_malla);
x_int = rho_int.*cos(theta_malla);
y_int = rho_int.*sin(theta_malla);
axis([-2, 2, -2, 2]); %AJUSTE DE LOS EJES
axis equal;
plot(x_ext, y_ext, 'b', 'LineWidth', 3); %SE REPRESENTA EL CILINDRO EXTERIOR
plot(x_int, y_int, 'b', 'LineWidth', 3); %SE REPRESENTA EL CILINDRO EXTERIOR
view(2) %VISTA EN 2D
hold off
2.1 Velocidad de las partículas del fluido y particularización de la ecuación
La velocidad de las partículas del fluido viene dada por la expresión [math]\vec{u}(\rho, \theta) = f(\rho)\vec{e}_{\theta}[/math], donde la presión p es constante. Sabiendo que (\vec{u},p) cumple la ecuacion de Navier-Stokes estacionaria se obtiene:
Siendo \mu el coeficiente de viscosidad del fluido.
Además, como la presión p es constante, su gradiente sera nulo: [math]\nabla p = 0[/math], y si aparte despreciamos el primer término de la fórmula anterior obtenemos que:
Así, se llega al Laplaciano de un campo vectorial que se desarrollara en la base cilíndrica de acuerdo a la siguiente expresión:
2.2 Cálculo del gradiente de una divergencia
Se desarrolla el primer termino, que es el gradiente de una divergencia en coordenadas cilindricas:
Si se particulariza al ejercicio se obtiene que:
Llegandose asi a la conclusión de que el primer término es igual al vector nulo: [math]\nabla(\nabla \cdot \vec{u}) = \vec{0}[/math]
2.3 Cálculo del rotacional del rotacional
Ahora se desarrolla el segundo termino, el rotacional del rotacional de [math]\vec{u}[/math], para ello se empieza en primer lugar con el primer rotacional:
Así, particularizando al ejercicio se obtiene:
Tras ello, se desarrolla el siguiente rotacional:
Por tanto, tras desarrollar todos los terminos se llega a que:
Ya, una vez conocidos todos los términos de la ecuación de Navier-Stokes, se puede resolver:
De lo que se deduce la siguiente ecuación diferencial:
Primero se despeja [math]f(\rho)[/math] y luego se simplifica llegando a la ecuación final:
2.5 Comprobación de la ecuación obtenida
A continuación se comprueba si la ecuación obtenida cumple la solución obtenida:
Para ello hay que derivar la expresión obtenida:
Tras haber realizado la derivada y sustituir en la ecuación, se comprueba que la ecuación cumple la solución obtenida:
2.6 Calculo de los parámetros a y b
Ahora queda determinar a y b para que las velocidades en la frontera del fluido coincidan con las de los cilindros. Para ello, se aplican las condiciones obtenidas sobre los parámetros obtenidos. Primero se particulariza en el cilindro interior, donde [math]\rho[/math]=1:
Segundo se particulariza en el cilindro exterior, donde [math]\rho[/math]=2:
Solucionando el sistema en función de las velocidades angulares se obtiene que:
Tras averiguar los valores de a y b, [math]\vec{u}[/math] queda así:
2.7 Condición de incompresibilidad
Para comprobar la condición de incompresibilidad, se hace la divergencia del [math]\vec{u}[/math] obtenido en coordenadas cilindricas y se comprueba que la divergencia del campo de velocidades es nulo:
3 Representación del campo de velocidades
Partiendo de las siguientes expresiones; [math]\vec{u}(\rho,\theta) = f(\rho)\,\vec{e}_\theta,\quad f(\rho) = a\rho + \frac{b}{\rho}[/math]; y los valores anteriormente calculados; [math]b = \frac{2}{3} \omega_e + \frac{4}{3} \omega_i[/math], [math]a = - \frac{2}{3} \omega_e - \frac{1}{3} \omega_i[/math], si suponemos que [math]|\vec{\omega}_e| = |\vec{\omega}_i| = 1 \text{ y } \mu = 1[/math], podemos representar graficamente nuestro campo de velocidades:
clear; clc;
Ri = 1; % radio interior
Re = 2; % radio exterior
omega_i = +1; % interior: antihorario (convención +)
omega_e = -1; % exterior: horario (convención -)
B = (omega_i - omega_e) * (Ri^2 * Re^2) / (Re^2 - Ri^2);
A = omega_i - B / (Ri^2);
Ntheta = 20;
Nr = 15;
thetas = linspace(0, 2*pi, Ntheta+1);
thetas(end) = []; % evita duplicar 0 y 2*pi
rhos = linspace(Ri, Re, Nr);
figure('Color',[1 1 1]); hold on;
scale_quiver = 0.7;
for th = thetas
% velocidad tangencial u_theta en cada rho
u_theta = A .* rhos + B ./ rhos;
% componentes cartesianas (tangenciales)
U = -u_theta .* sin(th);
V = u_theta .* cos(th);
% posiciones
X = rhos .* cos(th);
Y = rhos .* sin(th);
% vectores del campo
quiver(X, Y, U, V, scale_quiver, 'LineWidth', 1.1, ...
'MaxHeadSize', 1.6, 'Color','r'); % azul verdoso
end
tc = linspace(0,2*pi,500);
plot(Ri*cos(tc), Ri*sin(tc), 'k', 'LineWidth', 2.0);
plot(Re*cos(tc), Re*sin(tc), 'k', 'LineWidth', 2.0);
axis equal;
xlim([-2.3, 2.3]); ylim([-2.3, 2.3]);
xlabel('X (m)');
ylabel('Y (m)');
title('Gráfica del campo de velocidades');
grid on;
holf off;
4 Líneas de Corriente
4.1 Campo v ortogonal a u
[math]\vec{v}(\rho, \theta) = \left[ \left( \left(-\frac{2}{3}w_{e}-\frac{1}{3}w_{i}\right)\rho \right) + \left( \frac{2}{3}w_{e}+\frac{4}{3}w_{i} \right) \frac{1}{\rho} \right] \vec{e}_{\rho}[/math]
4.2 Irrotacionalidad de v
[math]\text{Campo vectorial: } \vec{v}(\rho, \theta) = f(\rho)\vec{e}_{\rho} \implies v_{\rho} = f(\rho), \ v_{\theta} = 0, \ v_{z} = 0[/math]
[math]\text{Fórmula del Rotacional en Coordenadas Cilíndricas:}[/math]
[math]\nabla \times \vec{v} = \frac{1}{\rho}\left[\frac{\partial v_{z}}{\partial \theta} - \frac{\partial (\rho v_{\theta})}{\partial z}\right]\vec{e}_{\rho} + \left[\frac{\partial v_{\rho}}{\partial z} - \frac{\partial v_{z}}{\partial \rho}\right]\vec{e}_{\theta} + \frac{1}{\rho}\left[\frac{\partial (\rho v_{\theta})}{\partial \rho} - \frac{\partial v_{\rho}}{\partial \theta}\right]\vec{e}_{z}[/math]
\text{Sustitución de Componentes y Derivadas:}
[math]* \text{Componente } \vec{e}_{\rho}: \frac{1}{\rho}\left[\frac{\partial (0)}{\partial \theta} - \frac{\partial (\rho \cdot 0)}{\partial z}\right] = 0[/math]
[math]* \text{Componente } \vec{e}_{\theta}:
\left[\frac{\partial f(\rho)}{\partial z} - \frac{\partial (0)}{\partial \rho}\right] = 0 - 0 = 0 \quad (\text{porque } f(\rho) \text{ no depende de } z)[/math]
[math]* \text{Componente } \vec{e}_{z}:
\ltmath\gt\frac{1}{\rho}\left[\frac{\partial (\rho \cdot 0)}{\partial \rho} - \frac{\partial f(\rho)}{\partial \theta}\right] = \frac{1}{\rho}[0 - 0] = 0 \quad (\text{porque } f(\rho) \text{ no depende de } \theta)[/math]
\text{Conclusión: El campo es irrotacional.}
[math]\nabla \times \vec{v} = \vec{0}[/math]
4.3 Función de corriente de u (𝜓)
La función de corriente de 𝑢⃗ es el potencia escalar de dicho campo, lo que significa: ∇𝜓 = v
[math]\nabla \psi = \frac{\partial \psi}{\partial \rho} \vec{e}_{\rho} + \frac{1}{\rho} \frac{\partial \psi}{\partial \theta} \vec{e}_{\theta} = \vec{v}[/math]
[math]\frac{\partial \psi}{\partial \rho} = \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \rho + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \frac{1}{\rho}[/math]
[math]\frac{1}{\rho} \frac{\partial \psi}{\partial \theta} = 0 \implies \frac{\partial \psi}{\partial \theta} = 0[/math]
Podemos ver que ψ depende únicamente de ρ, por lo que integramos la ecuación respecto a ρ
[math]\psi(\rho) = \int \left[ \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \rho + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \frac{1}{\rho} \right] d\rho[/math]
[math]\psi(\rho) = \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \int \rho \, d\rho + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \int \frac{1}{\rho} \, d\rho[/math]
[math]\psi(\rho) = \left(-\frac{2}{3}w_e - \frac{1}{3}w_i\right) \frac{\rho^2}{2} + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \ln(\rho) + C[/math]
El resultado sería:
[math]\psi(\rho, \theta) = \left(-\frac{1}{3}w_e - \frac{1}{6}w_i\right) \rho^2 + \left(\frac{2}{3}w_e + \frac{4}{3}w_i\right) \ln(\rho) + C[/math]
4.4 Comprobación de que son líneas de corriente
Como ya mencionamos anteriormente, para que v represente las líneas de correinte de u estos dos campos deben ser ortogonales entre sí.
[math]\vec{u} \cdot \nabla \psi = 0 \implies \vec{u} \cdot \vec{v} = 0[/math]
Componentes del campo tangencial [math]\vec{u}:[/math]
[math]\quad u_{\rho}=0, \quad u_{\theta}=f(\rho), \quad u_{z}=0[/math]
Componentes del radial ortogonal [math]\vec{v}:[/math] [math]\quad v_{\rho}=f(\rho), \quad v_{\theta}=0, \quad v_{z}=0[/math]
La multiplicación escalar sin sustituir los valores quedaría tal que asi:
[math]\vec{u} \cdot \vec{v} = u_{\rho}v_{\rho} + u_{\theta}v_{\theta} + u_{z}v_{z}[/math]
Procedemos a sustituir las componentes por las aportadas anteriormente:
[math]\vec{u} \cdot \vec{v} = (0) \cdot (f(\rho)) + (f(\rho)) \cdot (0) + (0) \cdot (0)[/math]
[math]\vec{u} \cdot \vec{v} = 0 + 0 + 0 = 0[/math]
Como podemos ver se cumple la ortogonalidad ya que el resultado es = 0
4.5 Líneas de corriente
%% Parámetros
Ri = 1; Re = 2;
we = -1; wi = +1;
% Coeficientes de la función f(rho) = A_f*rho + B_f/rho
A_f = (-2/3)*we - (1/3)*wi;
B_f = (2/3)*we + (4/3)*wi;
A_Psi = A_f / 2; % Coeficiente de R^2
B_Psi = B_f; % Coeficiente de log(R)
%% Mallado y Cálculo de Psi
L = 3; N = 300;
rho = linspace(0.001, L, N);
theta = linspace(0, 2*pi, N);
[R,T] = meshgrid(rho, theta);
X = R.*cos(T);
Y = R.*sin(T);
% Cálculo del potencial escalar
Psi = A_Psi*(R.^2) + B_Psi*log(R);
% Enmascaramiento de zonas fuera del anillo
mask = (R < Ri | R > Re);
Psi(mask) = NaN;
%% Representación Gráfica
figure;
hold on;
% Dibujar las líneas de corriente
contour(X, Y, Psi, 20, 'r', 'LineWidth', 1.5);
% Dibujar los cilindros (fronteras)
th = linspace(0, 2*pi, 400);
plot(Ri*cos(th), Ri*sin(th), 'k', 'LineWidth', 2);
plot(Re*cos(th), Re*sin(th), 'k', 'LineWidth', 2);
% Configuración del gráfico
axis equal;
axis([-L, L, -L, L]);
xlabel('x');
ylabel('y');
title('\Psi(\rho) = cte. Líneas de Corriente del campo U');
grid on;
hold off;
<pre>
}}
==Puntos de velocidad máxima del fluido==
Para enconntrar los puntos dentro de los cilindros donde la velocidad de nuestro fluido es máxima, debemos estudiar la primera derivada del campo de velocidades obtenido: <math>\vec{u}(\rho,\theta) = \left[ \left( -\frac{2}{3}\omega_e - \frac{1}{3}\omega_i \right)\rho + \frac{\frac{2}{3}\omega_e + \frac{4}{3}\omega_i}{\rho} \right]\vec{e}_\theta</math>
== Rotacional de u ==
=== Cálculo del Rotacional ===
Las componentes del campo <math>u</math> son: <math>u_\rho = 0, u_\theta = f(\rho), u_z = 0</math>
La fórmula del rotacional de un vector u es: <math>\nabla \times \mathbf{u} = \frac{1}{\rho} \left[ \frac{\partial u_z}{\partial \theta} - \frac{\partial (\rho u_\theta)}{\partial z} \right] \mathbf{e}_\rho + \left[ \frac{\partial u_\rho}{\partial z} - \frac{\partial u_z}{\partial \rho} \right] \mathbf{e}_\theta + \frac{1}{\rho} \left[ \frac{\partial (\rho u_\theta)}{\partial \rho} - \frac{\partial u_\rho}{\partial \theta} \right] \mathbf{e}_z</math>
Al sustituir las componentes en la ecuación del rotacional de la componente <math>\vec{e}_{\rho}</math> nos queda lo siguiente:
<math>\frac{1}{\rho} \left[ \frac{\partial (0)}{\partial \theta} - \frac{\partial (\rho f(\rho))}{\partial z} \right]</math>
En este caso el resultado sería = 0 y con la componente <math>\vec{e}_{\theta}</math> ocurre lo mismo.
Por último calculamos la componente <math>\vec{e}_{z}</math> : <math>
(\nabla \times \mathbf{u})_{\mathbf{e}_z} = \frac{1}{\rho} \left[ \frac{\partial (\rho f(\rho))}{\partial \rho} - \frac{\partial (0)}{\partial \theta} \right]
</math>
El segundo valor es 0 así que solo debemos calcular el primero con <math>
\rho f(\rho) = \rho \left[ \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) \rho + \left( \frac{2}{3} w_e + \frac{4}{3} w_i \right) \frac{1}{\rho} \right] = \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) \rho^2 + \left( \frac{2}{3} w_e + \frac{4}{3} w_i \right)
</math>
Calculamos la derivada de eso: <math>
\frac{d}{d\rho}(\rho f(\rho)) = \frac{d}{d\rho} \left[ \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) \rho^2 + \text{Constante} \right] = \left( -\frac{2}{3} w_e - \frac{1}{3} w_i \right) (2\rho) + 0 = \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) \rho
</math>
Ahora solo debemos sustituir en la componente <math>\vec{e}_{z}</math>
<math>
(\nabla \times \mathbf{u})_{\mathbf{e}_z} = \frac{1}{\rho} \left[ \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) \rho \right] = \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right)
</math>
Por lo tanto, el rotacional de u es <math>
\nabla \times \mathbf{u} = \left( -\frac{4}{3} w_e - \frac{2}{3} w_i \right) \mathbf{e}_z
</math>
=== Representación del rotacional ===
[[Archivo:Rotacional5.jpg|700px|miniaturadeimagen|right|Leyenda: Rotacional de u.]]
{{matlab|codigo=
<pre>
%% Parámetros y Cálculo de la Magnitud Escalar
Ri = 1; Re = 2;
we = -1; wi = +1;
% La magnitud del rotacional es un escalar constante:
Rot = abs((-4/3) * we - (2/3) * wi);
% Calcular el valor numérico
fprintf('La magnitud del rotacional es constante: %.4f\n', Rot);
%% Mallado
L = 3; N = 300;
rho = linspace(0.001, L, N);
theta = linspace(0, 2*pi, N);
[R,T] = meshgrid(rho, theta);
X = R.*cos(T);
Y = R.*sin(T);
% Debe estar en una matriz para la función pcolor
Rot = ones(size(R)) * MagRot_val;
% Enmascaramos para solo mostrar el anillo del rotacional
mask = (R < Ri | R > Re);
Rot(mask) = NaN;
%% Representación Gráfica
figure;
hold on;
grid on;
% Usamos pcolor para mapear la magnitud del campo
pcolor(X, Y, Rot);
shading flat; % Suaviza la cuadrícula
colormap('hot'); % Utiliza un mapa de colores calido
% Dibujar los bordes
th = linspace(0, 2*pi, 400);
plot(Ri*cos(th), Ri*sin(th), 'k', 'LineWidth', 2);
plot(Re*cos(th), Re*sin(th), 'k', 'LineWidth', 2);
% Configuración del gráfico
axis equal;
axis([-L, L, -L, L]);
xlabel('x');
ylabel('y');
title('Rotacional de u');
5 Gradiente de la temperatura
5.1 =Representación grafica
En la figura 8 se puede apreciar como las curvas de nivel que representan la temperatura son ortogonales al gradiente.
clear all % Limpiar la pantalla y espacio de trabajo
clc
r=1:0.1:2; % parámetros
th=0:0.1:2*pi;
[R,theta]=meshgrid(r,th);
temp=log10(1+R.^2).*(((cos(theta)).^2).*exp(-(R-(5/2)).^2)); % Expresión de la temperatura
[Dx,Dy]=gradient(temp,r,th); % Gradiente de la temperatura
figure; % Crea una figura y dibuja las curvas de nivel
quiver(R,theta,Dx,Dy,'color','r'); % Dibuja un campo de vectores sobre las curvas de nivel
xlabel('r'), ylabel('\theta'), colorbar;
title('curvas de nivel');
hold on
contour(R,theta,temp,'b');
legend('vectores', 'curvas de nivel'); % Leyenda para indicar el campod e vectores
hold off5.2 Calculo del gradiente
El gradiente de la temperatura en coordenadas cilíndricas queda definido por la ecuacuacion:
[math]\nabla T(\rho,\theta) = \frac{\partial T}{\partial \rho}\,\vec e_{\rho} + \left(\frac{1}{\rho}\frac{\partial T}{\partial \theta}\right)\vec e_{\theta} + \frac{\partial T}{\partial z}\,\vec e_{z}[/math]
Cada coordenada tendrá su derivada parcial respecto a ella misma,las cuales sumadas son el gradiente:
6 Caudal de una sección longitudinal
Para comenzar calcularemos el caudal que atraviesa los dos cilindros por una sección exterior del cilindro.
Para ello tendremos en cuenta que los dos cilindros tiene 2 metros de profundidad y la velocidad del fluido esta condicionada por la ecuacion: [math]\vec{u}(\rho,\theta) = \left(\frac{4}{3}\rho - \frac{4}{3}\rho\right)\,\vec e_{\theta} [/math]
Donde depende de la variable se obtienen estos datos:
[math]\quad u_{\rho}=0,\quad u_{\theta}=\frac{4}{3}\rho - \frac{4}{3}\rho,\quad u_{z}=0 [/math]
Para calcular el caudal que recorre la superficie del cilindro(Imagen 9) se usa esta formula la cual consiste en integrar la derivada de la superficie por la velocidad:
[math] \int_{S} \vec{u}\cdot d\vec{S} \;=\; \iint_{D} \vec{u}\bigl(\Phi(u,v)\bigr)\,\cdot\, \left(\Phi_{u}\times\Phi_{v}\right)\,du\,dv [/math]
si elk caudal sale negativo significa que hemos tomado la dirreccion contraria al vector theta.
6.1 Parametrización de la superficie
Para poder calcular el caudal es necesario parametrizas la superficie por la cual circula el caudal (Imagen 9)
Esta es una de las posibles parametrizaciones de la superficie S:
[math] \vec{r}(\theta, z) = \begin{pmatrix} R \cos\theta \\ R \sin\theta \\ z \end{pmatrix}, \qquad 0 \le \theta \le 2\pi,\quad 0 \le z \le H. [/math]
calculando con esta parametrizacion nos da que el flujo es negativo por lo tanto buscamos otra que nos lo de correctamente:
6.2 Calculo del Caudal
Una vez obtenida la parametrizacion se puede obtener el vector u en funcion de esta:
R_exterior = 2.0; % Radio exterior
altura = 2.0;
% Generar cilindro exterior
[Xe, Ye, Ze] = cylinder(R_exterior, 80);
Ze = Ze * altura;
% GRAFICA
figure;
hold on;
grid on;
axis equal;
color_superficie = [0.1 0.7 0.8];
% Superficie cilíndrica exterior
surf(Xe, Ye, Ze, ...
'FaceColor', color_superficie, ...
'FaceAlpha', 0.7, ...
'EdgeColor', 'k');
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Pared exterior del cilindro');
view(45, 25);
hold off;
