Temperaturas, ondas y campos vectoriales en un arco (Grupo 28)

Trabajo realizado por estudiantes
Título	Temperaturas, ondas y campos vectoriales en un arco. Grupo 28
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Tiago di Risio Diego Gonzalez Ramirez Lucas Escalante Morante Nicolás Bofarull Esteban Alba García Celdrán
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Nuestro proyecto trabaja con un campo vectorial de un sector anular. Esta es una curva plana comprendida en el plano X-Y, por lo que su valor de Z siempre va a ser nulo (Z=0). Por otra parte la ρ esta comprendida entre 1 y 2 (ρ ∈[1, 2]), y Theta oscila de 0 a π (θ ∈[0, π]), por lo que seria como la sección horizontal de medio donut, o una semicircunferencia truncada en el centro.

1 Póster

Se puede acceder al póster esquema del trabajo a través del siguiente enlace:

https://upm365-my.sharepoint.com/:b:/g/personal/alba_gceldran_alumnos_upm_es/ETANcBmyHXhNiDgICsb3QpsB_ME3k-oPZ5WYsfSolvbYJA?e=ngMgb0

2 Representación del mallado

2.1 Código y representación

Representación del mallado

clear; clc; close all;

% Definición de parámetros
% Radio entre 1 y 2
rho_min = 1; 
rho_max = 2;

% Ángulo: Según la Figura 3, es un semicírculo (0 a pi)
theta_min = 0; 
theta_max = pi;

% Paso de muestreo h = 1/10 (0.1)
h = 0.1;

%% 2. Generación de la Malla (De polares a Cartesianas)
% Creamos los vectores de radio y ángulo
rho_vec = rho_min:h:rho_max;
theta_vec = [theta_min:h:theta_max, theta_max];

% Creamos la rejilla en coordenadas polares
[R, Th] = meshgrid(rho_vec, theta_vec);

% Transformamos a Coordenadas Cartesianas para poder pintar
% x = rho * cos(theta)
% y = rho * sin(theta)
X = R .* cos(Th);
Y = R .* sin(Th);

%% 3. Visualización (Replicando Figura 3)
figure(1);
clf; % nuevo proceso
hold on;
axis equal; % Mantiene las proporciones entre ejes
grid off; % Quitamos la rejilla de fondo.

% A) Pintar la malla interna (Verde)
% Pintamos las líneas radiales 
plot(X, Y, 'g', 'LineWidth', 0.5); 
% Pintamos los arcos 
plot(X', Y', 'g', 'LineWidth', 0.5);

% B) Pintar el Borde Rojo
% El borde se compone de 4 partes:
% 1. Arco interior (rho = 1)
% 2. Arco exterior (rho = 2)
% 3. Cierre inferior izquierdo y derecho

% Definimos los bordes para pintarlos seguidos
borde_rho = [rho_vec(1)*ones(size(theta_vec)), ...      % Arco Interior
             rho_vec(end)*ones(size(theta_vec)), ...    % Arco Exterior
             rho_vec, ...                               % Línea base (theta=0)
             rho_vec];                                  % Línea base (theta=pi)

% Para pintar el borde continuo, lo hacemos paramétricamente:
% Lado 1: Arco grande (Exterior)
plot(rho_max * cos(theta_vec), rho_max * sin(theta_vec), 'r', 'LineWidth', 2);
% Lado 2: Arco pequeño (Interior)
plot(rho_min * cos(theta_vec), rho_min * sin(theta_vec), 'r', 'LineWidth', 2);
% Lado 3: Cierre plano izquierda (theta = pi)
plot([rho_min*cos(pi), rho_max*cos(pi)], [rho_min*sin(pi), rho_max*sin(pi)], 'r', 'LineWidth', 2);
% Lado 4: Cierre plano derecha (theta = 0)
plot([rho_min*cos(0), rho_max*cos(0)], [rho_min*sin(0), rho_max*sin(0)], 'r', 'LineWidth', 2);

%% 4. Ajustes de la representacion
title('Mallado del Arco I');
xlabel('x');
ylabel('y');
axis([-2.5 2.5 -0.5 2.5]); % Ajustar límites para ver bien todo
set(gcf, 'Color', 'w'); % Fondo blanco

3 Temperatura del sólido

La temperatura del sólido se supone conocida y viene dada por la función:

[math]T(x,y)=(x-y)^2 [/math]

3.1 Código

Representación de las temperaturas

En la representación de la temperatura del arco, se observan las distintas líneas de nivel de la función temperatura con distintos colores, siendo los mas oscuros y fríos los de las temperaturas mas bajas y los mas brillantes y cálidos los de las mas altas.

clear; clc; close all;

% Definición de parámetros
% Radio entre 1 y 2
rho_min = 1; 
rho_max = 2;

% Ángulo: Según la Figura 3, es un semicírculo (0 a pi)
theta_min = 0; 
theta_max = pi;

% Paso de muestreo h = 1/10 (0.1)
h = 0.1;

%% 2. Generación de la Malla (De polares a Cartesianas)
% Creamos los vectores de radio y ángulo
rho_vec = rho_min:h:rho_max;
theta_vec = [theta_min:h:theta_max, theta_max];

% Creamos la rejilla en coordenadas polares
[R, Th] = meshgrid(rho_vec, theta_vec);

% Transformamos a Coordenadas Cartesianas para poder pintar
% x = rho * cos(theta)
% y = rho * sin(theta)
X = R .* cos(Th);
Y = R .* sin(Th);

%% 3. Visualización
figure(1);
clf; % nuevo proceso
hold on;
axis equal; % Mantiene las proporciones entre ejes
grid off; % Quitamos la rejilla de fondo.

% A) Pintar la malla interna (Verde)
% Pintamos las líneas radiales 
plot(X, Y, 'g', 'LineWidth', 0.5); 
% Pintamos los arcos 
plot(X', Y', 'g', 'LineWidth', 0.5);

% B) Pintar el Borde Rojo
% El borde se compone de 4 partes:
% 1. Arco interior (rho = 1)
% 2. Arco exterior (rho = 2)
% 3. Cierre inferior izquierdo y derecho

% Definimos los bordes para pintarlos seguidos
borde_rho = [rho_vec(1)*ones(size(theta_vec)), ...      % Arco Interior
             rho_vec(end)*ones(size(theta_vec)), ...    % Arco Exterior
             rho_vec, ...                               % Línea base (theta=0)
             rho_vec];                                  % Línea base (theta=pi)

% Para pintar el borde continuo, lo hacemos paramétricamente:
% Lado 1: Arco grande (Exterior)
plot(rho_max * cos(theta_vec), rho_max * sin(theta_vec), 'r', 'LineWidth', 2);
% Lado 2: Arco pequeño (Interior)
plot(rho_min * cos(theta_vec), rho_min * sin(theta_vec), 'r', 'LineWidth', 2);
% Lado 3: Cierre plano izquierda (theta = pi)
plot([rho_min*cos(pi), rho_max*cos(pi)], [rho_min*sin(pi), rho_max*sin(pi)], 'r', 'LineWidth', 2);
% Lado 4: Cierre plano derecha (theta = 0)
plot([rho_min*cos(0), rho_max*cos(0)], [rho_min*sin(0), rho_max*sin(0)], 'r', 'LineWidth', 2);

%% 4. Ajustes de la representacion
title('Mallado del Arco I');
xlabel('x');
ylabel('y');
axis([-2.5 2.5 -0.5 2.5]); % Ajustar límites para ver bien todo
set(gcf, 'Color', 'w'); % Fondo blanco

%% 5. Campo de Temperaturas
% Definimos la función T = (x - y)^2
T = (X - Y).^2;

% Creamos una nueva figura para no borrar la del mallado limpio
figure(2);
clf; hold on; axis equal;

% Mapa de Calor
[C, h_cont] = contourf(X, Y, T, 20, 'LineStyle', 'none'); 

% B) Añadir la Barra de Color
colorbar; 
title('Distribución de Temperatura T(x,y) = (x-y)^2');
xlabel('x'); ylabel('y');

% C) Añadir el Borde Negro (Contorno del arco)
plot(rho_max * cos(theta_vec), rho_max * sin(theta_vec), 'k', 'LineWidth', 2);
plot(rho_min * cos(theta_vec), rho_min * sin(theta_vec), 'k', 'LineWidth', 2);
plot([rho_min*cos(pi), rho_max*cos(pi)], [rho_min*sin(pi), rho_max*sin(pi)], 'k', 'LineWidth', 2);
plot([rho_min*cos(0), rho_max*cos(0)], [rho_min*sin(0), rho_max*sin(0)], 'k', 'LineWidth', 2);

Nuestro trabajo explicaba que tenemos que seguir el mismo proceso que en el K, con la diferencia de que nos dan una ecuación de temperatura distinta. En el K también indica que existe un foco de calor en rho igual a 1. En nuestra ecuación de temperatura eso no se cumple ya que es la indicada en el punto 2. Esta fórmula explica que la temperatura aumenta cuando la diferencia absoluta de la x y la y incrementa exponencialmente elevada a dos, explicado de una manera mas simple, la temperatura crece exponencialmente según se aleja de la línea x=y, en esa línea la temperatura siempre será 0.

Para visualizar de manera mas sencilla la forma en la que crece la temperatura según se aleja de la línea X=Y, representamos la función [math]T(x,y)=(x-y)^2 [/math] en Geogebra 3D de esta forma, se aprecia perfectamente como la función temperatura es un cilindro parabólico a lo largo del eje X=Y y con vértice en el plano Z=0.

• 

  Visualización de la forma de cilindro parabólico de la función

• 

  Visualización de la función proyectando el eje Z

4 Gradiente de T

4.1 Definición de un gradiente

El gradiente ([math]\nabla T[/math]) se utiliza para describir la dirección y tasa de cambio de más rápida de un campo escalar. El vector indica la dirección en la que varía más rápidamente y su módulo (|[math]\nabla T[/math]|) indica la tasa en esa dirección. Para cacular el gradiente en coordenadas cartesianas, se utiliza la siguiente expresión:

[math]\nabla T = \frac{\partial T}{\partial x}\vec i+\frac{\partial T}{\partial y}\vec j+\frac{\partial T}{\partial z}\vec k[/math]

Teniendo en cuenta la función de temperatura dada([math]T(x,y)=(x-y)^2 [/math]), el gradiente será:

[math]\nabla T = 2(x-y)\vec i-2(x-y)\vec j[/math]

Después de representar el gradiente de la función T sobre las líneas isotermas de la misma, se puede observar como el propio gradiente es perpendicular a dichas líneas en cada punto de la función.

4.2 Código y Representación

Representación del gradiente de T sobre las líneas isotermas

%% GRADIENTE DE TEMPERATURA
clear; clc; close all;

% --- 1. GEOMETRÍA ---
rho_vec = 1:0.05:2;
theta_vec = [0:0.05:pi, pi]; 
[R, Th] = meshgrid(rho_vec, theta_vec);

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

% --- 2. CÁLCULO DE CAMPOS ---
% Temperatura T = (x - y)^2
T = (X - Y).^2;

% Gradiente (Derivadas parciales)
% dT/dx = 2*(x - y)
% dT/dy = -2*(x - y)
TX = 2 * (X - Y);
TY = -2 * (X - Y);

% --- 3. VISUALIZACIÓN ---
figure(10); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Gradiente de Temperatura');
xlabel('x'); ylabel('y');

% Mapa de Color (Temperatura)
[C, h] = contourf(X, Y, T, 30, 'LineStyle', 'none');
c = colorbar;
ylabel(c, 'Temperatura T(x,y)');
colormap(parula); % Mapa de color estándar

% Flechas del Gradiente 
paso = 4; 
idx_r = 1:paso:size(X,1);
idx_t = 1:paso:size(X,2);

X_q  = X(idx_r, idx_t);
Y_q  = Y(idx_r, idx_t);
TX_q = TX(idx_r, idx_t);
TY_q = TY(idx_r, idx_t);

% flechas
quiver(X_q, Y_q, TX_q, TY_q, 'k', 'LineWidth', 1, 'AutoScaleFactor', 1);

% Bordes para que quede bonito
plot_borde(rho_vec, theta_vec);

axis([-2.5 2.5 0 2.5]); 
grid off;

% --- Función Borde ---
function plot_borde(r_v, t_v)
    col = 'k'; ancho = 2;
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

5 Campo de vectores

Dado el siguiente campo:

[math]\vec{u}(\rho,\theta)=\frac{1}{5}(\rho - 1)\rho^{2}\sin\theta\,\vec{e}_{\theta}[/math]

En la figura se puede ver con flechas rojas las componentes del campo vectorial. Las únicas representadas son las tangenciales, en otras palabras la [math]\vec{e}_{\theta}[/math]. La componente normal ([math]\vec{e}_{\rho}[/math]), y la componente binormal ([math]\vec{e}_{z}[/math]), son las dos nulas, iguales a 0, por eso mismo no tienen ninguna representación. La normal tendría una dirección alejándose o acercándose del centro del circulo dependiendo si es positiva o negativa. Y la componente binormal si todo fuese positivo se saldría de la pantalla hacia nosotros, direccion vertical. Estas tres componentes siempre son positivas y tienen que cumplir la regla de la mano derecha, cuando hablamos de sus orientaciones.

5.1 Código

Representación campo vectorial U

clear; clc; close all;

% 1. Definir Geometría 
rho_vec = 1:0.1:2;                % Radio de 1 a 2
theta_vec = 0:0.1:pi;             % De 0 a pi (Semicírculo)
[R, Th] = meshgrid(rho_vec, theta_vec); % Malla en polares

% Coordenadas 
X = R .* cos(Th);
Y = R .* sin(Th);

% 2. Calcular el Campo Vectorial u 
% Fórmula: u = 1/5 * (rho-1) * rho^2 * sin(theta) * e_theta
U_rho   = zeros(size(R));  % No hay componentes normales ni binormales
U_theta = (1/5) * (R - 1) .* (R.^2) .* sin(Th);

% Transformar vectores a Cartesianas 
UX = U_rho .* cos(Th) - U_theta .* sin(Th);
UY = U_rho .* sin(Th) + U_theta .* cos(Th);

% 3. Optimización visual 
paso = 2; % Pintar solo 1 de cada 2 flechas para que se vean nítidas
idx_r = 1:paso:size(X,1);
idx_t = 1:paso:size(X,2);

X_q  = X(idx_r, idx_t);
Y_q  = Y(idx_r, idx_t);
UX_q = UX(idx_r, idx_t);
UY_q = UY(idx_r, idx_t);

% 4. Pintar la Figura
figure(6); clf; hold on; axis equal;
set(gcf, 'Color', 'w'); % Fondo blanco
title('Campo Vectorial U');
xlabel('x'); ylabel('y');

% Pintar contorno del arco (Referencia visual)
borde_R = [1, 2, 2, 1, 1]; % Radios para dibujar el marco
borde_T = [0, 0, pi, pi, 0]; % Ángulos para dibujar el marco
% (Nota: pinto líneas simples de referencia)
plot(2*cos(0:0.01:pi), 2*sin(0:0.01:pi), 'k', 'LineWidth', 1.5); % Arco ext
plot(1*cos(0:0.01:pi), 1*sin(0:0.01:pi), 'k', 'LineWidth', 1.5); % Arco int
line([-2 -1], [0 0], 'Color', 'k', 'LineWidth', 1.5); % Cierre izq
line([1 2], [0 0], 'Color', 'k', 'LineWidth', 1.5);   % Cierre der

% Pintar las flechas 
quiver(X_q, Y_q, UX_q, UY_q, 'r', 'LineWidth', 1.5, 'MaxHeadSize', 0.5);

axis([-2.5 2.5 0 2.5]); % Ajustar zoom
grid off;

6 Dibujo del sólido antes y después del desplazamiento

Este apartado se dedica a la visualización geométrica del problema, mostrando cómo el campo vectorial de desplazamiento [math]\vec{u}[/math] transforma la forma original del sólido. La meta es comparar directamente el estado no deformado con el estado deformado. Se puede observar tanto una expansión que incrementa cuanto mas cerca esta del eje de simetría de nuestra figura, y también como se desplaza en sentido antihorario, cuyo desplazamiento incrementa según el ρ aumenta.

6.1 codigo

Inicial

Final

Comparación

%% Visualización de Deformación (Azul vs Rojo)
clear; clc; close all;
% --- 1. DATOS Y CÁLCULOS  ---
rho_vec = 1:0.1:2;

% EL CAMBIO ESTÁ AQUÍ:

theta_vec = [0:0.1:pi, pi]; 

[R, Th] = meshgrid(rho_vec, theta_vec);

% Posición Inicial
X_ini = R .* cos(Th);
Y_ini = R .* sin(Th);

% Desplazamiento u (Trabajo M)
u_rho   = zeros(size(R));
u_theta = (1/5) * (R - 1) .* (R.^2) .* sin(Th);

UX = u_rho .* cos(Th) - u_theta .* sin(Th);
UY = u_rho .* sin(Th) + u_theta .* cos(Th);

% Posición Final
X_fin = X_ini + UX;
Y_fin = Y_ini + UY;

% ---  GENERACIÓN DE LAS GRÁFICAS ---

% GRÁFICA 1: Posición Inicial
figure(1); clf; hold on; axis equal;
set(gcf, 'Color', 'w'); title('1. Posición Inicial (Sin deformar)');
xlabel('x'); ylabel('y');
plot(X_ini, Y_ini, 'b', 'LineWidth', 1);
plot(X_ini', Y_ini', 'b', 'LineWidth', 1);
plot_borde(rho_vec, theta_vec, 'k', 2); 
axis([-2.5 2.5 0 2.5]); grid on;

% GRÁFICA 2: Posición Final
figure(2); clf; hold on; axis equal;
set(gcf, 'Color', 'w'); title('2. Posición Final (Deformada)');
xlabel('x'); ylabel('y');
plot(X_fin, Y_fin, 'r', 'LineWidth', 1);
plot(X_fin', Y_fin', 'r', 'LineWidth', 1);
axis([-2.5 2.5 0 2.5]); grid on;

% GRÁFICA 3: Superposición (AZUL vs ROJO)
figure(3); clf; hold on; axis equal;
set(gcf, 'Color', 'w'); title('3. Comparativa: Inicial vs Final');
xlabel('x'); ylabel('y');

% A) Inicial: AZUL
plot(X_ini, Y_ini, 'b', 'LineWidth', 1);
plot(X_ini', Y_ini', 'b', 'LineWidth', 1);

% B) Final: ROJO
plot(X_fin, Y_fin, 'r', 'LineWidth', 1.2);
plot(X_fin', Y_fin', 'r', 'LineWidth', 1.2);


% --- Función para bordes ---
function plot_borde(r_v, t_v, col, ancho)
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

7 Divergencia

7.1 Definición de la divergencia

La divergencia de un campo vectorial ([math]\nabla \cdot \vec{u}[/math]) en un punto dado es una medida de la tasa a la que el flujo del campo se está expandiendo (saliendo) o contrayendo (entrando) en ese punto.

Es un valor escalar que te dice qué tan fuerte es una fuente o un sumidero de flujo en ese lugar. Para calcular la divergencia en coordenadas cilíndricas se utiliza la siguiente fórmula:

[math] \nabla \cdot \vec{U} = \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} (\rho U_{\rho}) + \frac{\partial U_{\theta}}{\partial \theta} + \frac{\partial}{\partial z} (\rho U_{z}) \right] [/math]

Reemplazando los valores del campo en las posiciones de U, obtenemos la siguiente expresión:

[math] \nabla \cdot \vec{U} = \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho} (0) + \frac{\partial}{\partial \theta} \left( \frac{1}{5} (\rho - 1)\rho^2 \sin\theta \right) + \frac{\partial}{\partial z} (0) \right] [/math]

El resultado final de la divergencia es el siguiente:

[math] \nabla \cdot \vec{U} = \frac{1}{5} (\rho - 1)\rho \cos\theta [/math]

7.2 Código y representación

Mapa de color de la divergencia

%% DIVERGENCIA 

clear; clc; close all;

% 1. Geometría
rho_vec = 1:0.05:2;          
theta_vec = [0:0.05:pi, pi]; 
[R, Th] = meshgrid(rho_vec, theta_vec);

% Paso a cartesianas solo para pintar (X, Y)
X = R .* cos(Th);
Y = R .* sin(Th);

% 2. Cálculo de la Divergencia
% Fórmula: (1/5) * (rho^2 - rho) * cos(theta)
Div = (1/5) * (R.^2 - R) .* cos(Th);

% 3. Visualización
figure(7); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Divergencia: Expansión y Compresión');
xlabel('x'); ylabel('y');

% mapa de colores
[C, h] = contourf(X, Y, Div, 30, 'LineStyle', 'none');

% Barra de color
cb = colorbar;
ylabel(cb, 'Cambio de Volumen');

% borde negro
plot_borde(rho_vec, theta_vec, 'k', 2);

% Definimos un mapa de colores "Divergente" (Rojo-Azul)
%Azul para compresión, Rojo para expansión
colormap(jet); 

axis([-2.5 2.5 0 2.5]); 
grid off;

% --- Función Borde  ---
function plot_borde(r_v, t_v, col, ancho)
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

Como la divergencia depende del coseno de Theta, este cambia de signo en π/2. Por este motivo en la parte derecha del grafico, la divergencia es positiva, experimentando así un aumento de volumen y en la parte izquierda, la divergencia toma valores negativos por lo que el volumen se contrae. Finalmente en la línea entorno a π/2 la divergencia es cercana a 0 por lo que prácticamente no hay cambios en el volumen.

8 Rotacional

El rotacional de un campo vectorial [math]|\nabla \times \vec{u}|[/math] es una operación que mide la tendencia de un campo a girar. Visualmente, puedes imaginar el rotacional introduciendo una pequeña rueda de paletas en el campo. Si el rotacional es distinto de cero [math]|\nabla \times \vec{u}|≠ 0[/math], la rueda girará, indicando vorticidad (rotación). Si el rotacional es cero [math]|\nabla \times \vec{u}|=0[/math], la rueda no girará. El campo se llama irrotacional o conservativo.

La fórmula resulta en un nuevo vector con componentes en las direcciones: [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]

Expandiendo el determinante, obtenemos las tres componentes: [math] \nabla \times \vec{U} = \left( \frac{1}{\rho}\frac{\partial U_{z}}{\partial \theta} - \frac{\partial U_{\theta}}{\partial z} \right)\vec{e}_{\rho} \;+\; \left( \frac{\partial U_{\rho}}{\partial z} - \frac{\partial U_{z}}{\partial \rho} \right)\vec{e}_{\theta} \;+\; \frac{1}{\rho} \left[ \frac{\partial}{\partial \rho}(\rho U_{\theta}) - \frac{\partial U_{\rho}}{\partial \theta} \right] \vec{e}_{z} [/math]

En el caso de nuestro campo, el rotacional es igual a la siguiente expresión:

[math] \nabla \times \vec{U} = \frac{1}{5} \sin(\theta) (4\rho^2 - 3\rho) \, \vec{e}_{z} [/math]

8.1 Código y representación

Mapa de color del Rotacional

%% ROTACIONAL
% Fórmula derivada analíticamente en cilíndricas:
% Rot_z = (1/rho) * d(rho*u_theta)/drho
% Resultado: (1/5) * (4*rho^2 - 3*rho) * sin(theta)

clear; clc; close all;

% 1. Definición de Geometría
rho_vec = 1:0.05:2;          
theta_vec = [0:0.05:pi, pi]; 
[R, Th] = meshgrid(rho_vec, theta_vec);

% Paso a cartesianas
X = R .* cos(Th);
Y = R .* sin(Th);

% 2. Cálculo del Rotacional (Magnitud en eje Z)
Rot = (1/5) * (4*(R.^2) - 3*R) .* sin(Th);

% 3. Visualización
figure(7); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Magnitud del Rotacional');
xlabel('x'); ylabel('y');

%mapa de calor
[C, h] = contourf(X, Y, Rot, 30, 'LineStyle', 'none');

% Barra de color
cb = colorbar;
ylabel(cb, 'Intensidad de Giro');

% Borde negro 
plot_borde(rho_vec, theta_vec, 'k', 2);

% Mapa de colores
colormap(jet); 

axis([-2.5 2.5 0 2.5]); 
grid off;

% --- Función Borde ---
function plot_borde(r_v, t_v, col, ancho)
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

El campo gira más intensamente donde la función [math]\sin\theta[/math] es máxima (en el centro) y donde el radio [math]\rho[/math] es máximo (en el borde exterior), debido a que la velocidad tangencial aumenta desproporcionadamente con la distancia.

9 Tensiones normales

El cálculo de las tensiones se basa en la Ley de Hooke para un medio elástico lineal e isótropo, que define el tensor de tensiones [math]\mathbf{\sigma}[/math] a partir del tensor de deformaciones [math]\mathbf{\epsilon}[/math] y el cambio de volumen ([math]\nabla \cdot \vec{u}[/math]): [math]\mathbf{\sigma} = \lambda (\nabla \cdot \vec{u}) \mathbf{I} + 2\mu \mathbf{\epsilon}[/math]

Donde [math]\mathbf{\lambda}[/math] (el coeficiente relacionado con la resistencia a la dilatación volumétrica) y [math]\mathbf{\mu}[/math] (el módulo de cizalladura o resistencia al corte) son los Coeficientes de Lamé. Para este análisis, se toma el caso simplificado donde [math]\mathbf{\lambda = 1}[/math] y [math]\mathbf{\mu = 1}[/math]. Simplificación de la Ley de Hooke Al sustituir [math]\lambda=1[/math] y [math]\mu=1[/math] en la fórmula general para las tensiones normales ([math]\sigma[/math]), esta se simplifica a:

[math]\mathbf{\sigma = (\nabla \cdot \vec{u}) + 2\epsilon}[/math]

En nuestro caso, dado que el campo deformación ([math]\vec{u}(\rho,\theta)=\frac{1}{5}(\rho - 1)\rho^{2}\sin\theta\,\vec{e}_{\theta}[/math]) solo tiene componente en [math]{e}_{\theta}[/math] , al resolver la ecuación de la Ley de Hooke salen los siguientes resultados: Tensión Normal Radial: [math]\mathbf{\sigma_{\rho\rho}} = \underbrace{\left[ \frac{1}{5} (\rho^2 - \rho) \cos\theta \right]}_{\nabla \cdot \vec{u}} + 2 \underbrace{\left[ 0 \right]}_{\epsilon_{\rho\rho}}[/math] [math]\mathbf{\sigma_{\rho\rho} = \frac{1}{5} (\rho^2 - \rho) \cos\theta}[/math]

Tensión Normal Tangencial: [math]\sigma_{\theta\theta} = \underbrace{\left[ \frac{1}{5} (\rho^2 - \rho) \cos\theta \right]}_{\nabla \cdot \vec{u}} + 2 \underbrace{\left[ \frac{1}{5} (\rho^2 - \rho) \cos\theta \right]}_{\epsilon_{\theta\theta}}[/math]

Donde el resultado final es: [math]\mathbf{\sigma_{\theta\theta} = \frac{3}{5} \left( 1 - \frac{1}{\rho} \right) \cos\theta}[/math]

Teniendo en cuenta las fórmulas indicadas donde hay que reemplazar por los ejes( [math]\vec{e}_{\rho}[/math] y [math]\frac{1}{\rho}\vec{e}_{\theta}[/math]), los resultados finales de las tensiones normales son los siguientes:

[math]\mathbf{\sigma_{\rho\rho}} = \vec{e}_{\rho} \cdot \mathbf{\sigma} \cdot \vec{e}_{\rho} = (\nabla \cdot \vec{u}) + 2\epsilon_{\rho\rho}[/math]

[math]\mathbf{\sigma_{\rho\rho}} = \frac{1}{5} (\rho^2 - \rho) \cos\theta[/math]

[math]\mathbf{\sigma_{\theta\theta}} = \vec{e}_{\theta} \cdot \mathbf{\sigma} \cdot \vec{e}_{\theta} = (\nabla \cdot \vec{u}) + 2\epsilon_{\theta\theta}[/math]

[math]\mathbf{\sigma_{\theta\theta}} = \frac{3}{5} (\rho^2 - \rho) \cos\theta[/math]

Vistos los resultados, la tensión normal tangencial es 3 veces mayor que la tensión normal radial.

9.1 Código y Representación

Representación Tensión Normal Tangencial

Representación Tensión Normal Radial

%% TENSIONES NORMALES
clear; clc; close all;

% --- 1. GEOMETRÍA Y DATOS ---
rho_vec = 1:0.05:2;
theta_vec = [0:0.05:pi, pi]; 
[R, Th] = meshgrid(rho_vec, theta_vec);

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

% --- 2. CÁLCULO DE DEFORMACIONES (STRAIN) ---
% u_theta = 1/5 * (rho^3 - rho^2) * sin(theta)
% Epsilon_theta (Deformación angular)
% Fórmula: (1/rho) * du_theta/dtheta + u_rho/rho
E_theta = (1/5) * (R.^2 - R) .* cos(Th);

% Epsilon_rho (Deformación radial)
% Como no hay movimiento radial (u_rho=0), la deformación es 0.
E_rho = zeros(size(R)); 

% Divergencia
Div = E_rho + E_theta;

% --- 3. CÁLCULO DE TENSIONES (STRESS) ---
% Coeficientes
lambda = 1; 
mu = 1;

% Ley de Hooke en Polares:
Sigma_rr = lambda * Div + 2 * mu * E_rho;   % Tensión Radial
Sigma_tt = lambda * Div + 2 * mu * E_theta; % Tensión Tangencial

% --- 4. VISUALIZACIÓN ---

% FIGURA 8: Tensión Radial
figure(8); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Tensión Normal Radial');
xlabel('x'); ylabel('y');
[C, h] = contourf(X, Y, Sigma_rr, 30, 'LineStyle', 'none');
c = colorbar; ylabel(c, 'Pascales (Pa)');
colormap(gca, winter); % Azul/Verde
plot_borde(rho_vec, theta_vec);
axis([-2.5 2.5 0 2.5]); grid off;

% FIGURA 9: Tensión Tangencial (Colores Cálidos)
figure(9); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Tensión Normal Tangencial');
xlabel('x'); ylabel('y');
[C, h] = contourf(X, Y, Sigma_tt, 30, 'LineStyle', 'none');
c = colorbar; ylabel(c, 'Pascales (Pa)');
colormap(gca, autumn); % Rojo/Amarillo
plot_borde(rho_vec, theta_vec);
axis([-2.5 2.5 0 2.5]); grid off;

% Función auxiliar para bordes
function plot_borde(r_v, t_v)
    col = 'k'; ancho = 1.5;
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

En este apartado hemos analizado las fuerzas de estiramiento y compresión dentro del arco. Lo más llamativo es que el material sufre tres veces más a lo largo de la curva (tensión tangencial) que a lo ancho (tensión radial). Esto nos dice que, si la pieza rompiera por exceso de carga, la grieta aparecería perpendicularmente a la dirección del arco.

Además, hemos comprobado algo curioso: aunque el arco no se hace más grueso ni más fino (no hay deformación radial), sí que existe presión en esa dirección. Esto se debe a que el material es elástico y está conectado: al estirarlo mucho en una dirección, se generan fuerzas internas en las otras.

10 Tensiones tangenciales respecto al plano ortogonal [math]\vec{e}_{\rho}[/math]

En este apartado calculamos las tensiones tangenciales que actúan sobre el plano ortogonal al vector radial [math]\vec{e}_{\rho}[/math]. Siguiendo las instrucciones generales, la definición vectorial de esta tensión tangencial (vector de tracción menos su proyección normal) es: [math]\mathbf{\tau} = \mathbf{\sigma} \cdot \vec{e}_{\rho} - (\vec{e}_{\rho} \cdot \mathbf{\sigma} \cdot \vec{e}_{\rho})\vec{e}_{\rho}[/math]

Al simplificar esta expresión en coordenadas polares, la componente normal se anula, y nos queda la magnitud de la tensión cortante o de cizalladura: [math]|\mathbf{\tau}| = |\mathbf{\sigma_{\rho\theta}}|[/math]

Aplicando la Ley de Hooke para la cizalladura [math]\mathbf{\sigma_{\rho\theta} = 2 \epsilon_{\rho\theta}}[/math] con [math]\mu=1[/math]

Sustituyendo la deformación angular calculada previamente, la expresión final a representar es:

[math]|\mathbf{\tau}| = \left| \frac{1}{5} (2\rho^2 - \rho) \sin\theta \right|[/math]

Se buscará en la gráfica el punto donde esta magnitud es máxima, indicando dónde el material sufre mayor riesgo de desgarro.

10.1 Código y Representación

Representación Tensión Tangencial de Corte

%% APARTADO 9: TENSIONES TANGENCIALES 
clear; clc; close all;

% --- 1. GEOMETRÍA ---
rho_vec = 1:0.01:2;         
theta_vec = [0:0.01:pi, pi]; 
[R, Th] = meshgrid(rho_vec, theta_vec);

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

% --- 2. CÁLCULO DE LA TENSIÓN (Tau) ---
% Fórmula derivada: (1/5) * (2*rho^2 - rho) * sin(theta)
Tau = (1/5) * (2*R.^2 - R) .* sin(Th);

% Calculamos el valor absoluto
Tau_Mag = abs(Tau); 

% --- 3. VISUALIZACIÓN ---
figure(9); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Apartado 9: Tensión Tangencial de Corte (\tau_{\rho\theta})');
xlabel('x'); ylabel('y');

% Mapa de calor
[C, h] = contourf(X, Y, Tau_Mag, 50, 'LineStyle', 'none');
c = colorbar;
ylabel(c, 'Esfuerzo de Corte (Pa)');

% Usamos 
colormap(gca); 

% Bordes
plot_borde(rho_vec, theta_vec, 'k', 2);


% Buscamos el máximo
max_val = max(Tau_Mag(:));
[fil, col] = find(Tau_Mag == max_val);
x_max = X(fil, col);
y_max = Y(fil, col);

% Marcamos el punto máximo 
plot(x_max, y_max, 'wx', 'LineWidth', 2, 'MarkerSize', 10);
text(x_max, y_max+0.2, ' Máx', 'Color', 'k', 'FontWeight', 'bold');

axis([-2.5 2.5 0 2.5]); 
grid off;

% Función auxiliar borde
function plot_borde(r_v, t_v, col, ancho)
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

Aquí hemos calculado la fuerza que intenta 'rasgar' el material o hacer que sus capas deslicen unas sobre otras. Hemos localizado el punto más peligroso en la parte superior central del arco y en el borde exterior ([math]\rho=2[/math]).Este resultado encaja perfectamente con lo que vimos en la gráfica de la deformación: esa zona es justo donde la rejilla original se retuerce más.

11 Tensiones tangenciales respecto al plano ortogonal [math]\frac{1}{\rho} \vec{e}_{\theta}[/math]

La fórmula es la siguiente: [math]\mathbf{\tau} = \mathbf{\sigma} \cdot \left(\frac{1}{\theta}\vec{e}_{\theta}\right) - \left[ \left(\frac{1}{\rho}\vec{e}_{\theta}\right) \cdot \mathbf{\sigma} \cdot \left(\frac{1}{\rho}\vec{e}_{\theta}\right) \right] \left(\frac{1}{\rho}\vec{e}_{\theta}\right)[/math]

La fórmula utilizada será la siguiente: [math]|\mathbf{\tau}| = \left| \mathbf{\sigma_{\rho\theta}} \right|[/math]. Al reemplazar y aplicar la Ley de Hooke que se aplico en el apartado anterior obtenemos esta expresión:

[math]|\mathbf{\tau}| = \left| \left[ \frac{1}{5} (2\rho^2 - \rho) \sin\theta \right] \right|[/math]

11.1 Código y Representación

Representación Tensiones tangenciales respecto al plano ortogonal [math]\frac{1}{\rho} \vec{e}_{\theta}[/math]

%% APARTADO 10: TENSIONES TANGENCIALES (Plano 1/rho * e_theta)

clear; clc; close all;

% --- 1. GEOMETRÍA ---
rho_vec = 1:0.01:2;
theta_vec = [0:0.01:pi, pi]; 
[R, Th] = meshgrid(rho_vec, theta_vec);

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

% --- 2. CÁLCULO DE LA TENSIÓN ---
% Por simetría de tensiones (Tau_theta_rho = Tau_rho_theta)
% Usamos la misma fórmula derivada analíticamente en el Ap. 9
Tau = (1/5) * (2*R.^2 - R) .* sin(Th);

% Magnitud
Tau_Mag = abs(Tau); 

% --- 3. VISUALIZACIÓN ---
figure(10); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Apartado 10: Tensión Tangencial \tau_{\theta\rho}');
xlabel('x'); ylabel('y');

% Mapa de calor 
[C, h] = contourf(X, Y, Tau_Mag, 50, 'LineStyle', 'none');
c = colorbar;
ylabel(c, 'Esfuerzo de Corte (Pa)');
colormap(gca); 

% Bordes
plot_borde(rho_vec, theta_vec);

% Marcamos el máximo
max_val = max(Tau_Mag(:));
[fil, col] = find(Tau_Mag == max_val);

plot(X(fil, col), Y(fil, col), 'wx', 'LineWidth', 2, 'MarkerSize', 10);
text(X(fil, col), Y(fil, col)+0.2, ' Máx', 'Color', 'k', 'FontWeight', 'bold');

axis([-2.5 2.5 0 2.5]); 
grid off;

% --- Función Auxiliar ---
function plot_borde(r_v, t_v)
    col = 'k'; ancho = 1.5;
    plot(r_v(end)*cos(t_v), r_v(end)*sin(t_v), col, 'LineWidth', ancho);
    plot(r_v(1)*cos(t_v), r_v(1)*sin(t_v), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(1)), [r_v(1) r_v(end)]*sin(t_v(1)), col, 'LineWidth', ancho);
    plot([r_v(1) r_v(end)]*cos(t_v(end)), [r_v(1) r_v(end)]*sin(t_v(end)), col, 'LineWidth', ancho);
end

En este apartado hemos comprobado que la tensión de corte nos da exactamente lo mismo que en el apartado anterior. Esto tiene todo el sentido físico: para que el trozo de material se mantenga quieto y no se ponga a girar sobre sí mismo, las fuerzas que lo empujan por un lado tienen que ser iguales a las que lo empujan por el otro.

Básicamente, confirmamos que el punto más crítico donde el material intenta 'rasgarse' es la parte de arriba del arco y el borde exterior, que coincide justo con la zona donde vemos que los cuadrados de la malla se deforman más y se vuelven rombos.

12 Calculo de la masa

Procederemos con el calculo de la masa dada la función de la densidad

La masa en coordenadas polares se haya con la integral: [math] M = \int_{\theta_{\min}}^{\theta_{\max}} \int_{\rho_{\min}}^{\rho_{\max}} d(\rho,\theta)\, \rho \, d\rho \, d\theta [/math]

Donde:

La densidad está dada por: [math]d(\rho,\theta) = 1 + e^{\rho^{2}} \cos\theta[/math]

Y ρdρ es el jacobiano, que se podría calcular a mano también.

El dominio es un arco definido por:

[math]1 \le \rho \le 2, \qquad 0 \le \theta \le \pi[/math]

Por lo que serán los limites inferiores y superiores de la integral.

Aplicando la formula a nuestra densidad y los limites de integración:

[math] M = \int_{0}^{\pi} \int_{1}^{2} \left( 1 + e^{\rho^{2} \cos\theta} \right)\, \rho \, d\rho \, d\theta [/math]

Esta integral no se puede resolver a mano, así que usaremos aproximación numérica con la ayuda del programa Matlab.

% CÁLCULO MASA 
clear; clc;

% Función de densidad
% Multiplicamos por 'rho'(jacobiano)
fun = @(rho, theta) (1 + exp(rho.^2 .* cos(theta))) .* rho;

% Límites
rho_min = 1;  rho_max = 2;
theta_min = 0; theta_max = pi;

% Cálculo numérico de alta precisión
Masa = integral2(fun, rho_min, rho_max, theta_min, theta_max);

fprintf('La MASA TOTAL es: %.4f\n', Masa);

El resultado de la masa es: [math] M \approx 4.7124 + 19.9313 \approx \mathbf{24.6437} [/math]

13 Interpretación del trabajo

13.1 Motivación del Modelo

Interpretamos el dominio del arco como la sección transversal del revestimiento de hormigón (bóveda) de un túnel subterráneo. En ingeniería sísmica, una de las deformaciones más críticas que sufren los túneles es la distorsión por cortante , provocada por el paso de ondas sísmicas transversales (Ondas S) a través del terreno circundante.

13.2 Interpretación de los Resultados

El campo de desplazamientos impuesto en este trabajo, [math]\vec{u}(\rho,\theta)[/math], caracterizado por un movimiento puramente angular sin componente radial, simula eficazmente la respuesta cinemática del túnel ante este tipo de solicitación del terreno:

• Deformación (Apartados 4-7): La distorsión angular de la malla observada en las gráficas reproduce el comportamiento del anillo de hormigón al ser "arrastrado" diferencialmente por el suelo, perdiendo su ortogonalidad original.

• Tensiones Normales (Apartado 8): El hecho de que la tensión tangencial ([math]\sigma_{\theta\theta}[/math]) sea tres veces superior a la radial indica que el fallo principal del revestimiento se produciría por tracción/compresión excesiva en la fibra del arco, pudiendo generar grietas longitudinales a lo largo del túnel.

• Cizalladura (Apartados 9-10): Los máximos de tensión cortante ([math]\tau[/math]) localizados en la clave ([math]\theta=90^{\circ}[/math]) y en los laterales, nos alertan de puntos críticos donde el revestimiento del túnel podría sufrir deslizamientos relativos si no se disponen los conectores adecuados.

13.3 Conclusión de Diseño

El análisis sugiere que, para resistir estas cargas, el diseño del túnel no debe basarse solo en la resistencia a la compresión radial (presión de tierra), sino que requiere una armadura de acero robusta para soportar los esfuerzos de flexión y cortante inducidos por la deformación sísmica asimétrica.