Ondas en un arco y cálculo vectorial (66)

Ondas en un arco y cálculo vectorial

Trabajo realizado por estudiantes
Título	Ondas en un arco y cálculo vectorial. Grupo 66
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Nombres Jorge Lianes , Paula Calderón , Marina Morillo , Marta García-Inés y Sofía Romero
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Trabajo N: Arco II

Introducción

1 Dibujar un mallado que represente los puntos interiores del sólido

Mallado que representa los puntos interiores del sólido:

Se construye una rejilla rectangular con paso [math]h = \frac{1}{10}[/math] en [math]x[/math] e [math]y[/math]. Para cada línea vertical [math]x = x_i[/math] se calculan los límites interiores del semianillo como [math]y \in \big[\max\{0,\sqrt{1 - x_i^2}\}, \sqrt{4 - x_i^2}\big][/math], y se dibuja únicamente ese tramo. De forma análoga, para cada horizontal [math]y = y_j[/math] se dibujan los segmentos [math]x \in \big[\sqrt{1 - y_j^2}, \sqrt{4 - y_j^2}\big][/math] y su simétrico negativo, siempre que existan.

El dominio del semianillo está definido por la condición:

[math]1 \le \sqrt{x^2 + y^2} \le 2 \quad \text{y} \quad y \ge 0[/math]

Así el mallado queda aislado y representa exclusivamente los puntos interiores del sólido.

Código MATLAB y representación

Figura 1.1: Mallado puntos interiores

clc; clear; close all;
% Definición del dominio del semianillo
R_int = 1;      % radio interior
R_ext = 2;      % radio exterior
h = 0.1;        % paso de muestreo en x e y
% Rejilla cartesiana en el semiplano superior
x = -R_ext:h:R_ext;
y = 0:h:R_ext;
% Iniciar figura
figure; hold on; axis equal; grid on;
% Dibujar líneas verticales recortadas al semianillo
for i = 1:length(x)
    xi = x(i);
    if abs(xi) > R_ext, continue; end
    y_min_sq = R_int^2 - xi^2;
    y_max_sq = R_ext^2 - xi^2;
    if y_max_sq <= 0, continue; end
    y_min = max(0, sqrt(max(0, y_min_sq)));
    y_max = sqrt(y_max_sq);
    if y_max > y_min
        plot([xi xi], [y_min y_max], 'c');
    end
end
% Dibujar líneas horizontales recortadas al semianillo
for j = 1:length(y)
    yj = y(j);
    if yj > R_ext, continue; end
    x_min_sq = R_int^2 - yj^2;
    x_max_sq = R_ext^2 - yj^2;
    if x_max_sq <= 0, continue; end
    x_min = sqrt(max(0, x_min_sq));
    x_max = sqrt(x_max_sq);
    if x_max > x_min
        plot([-x_max -x_min], [yj yj], 'c');
        plot([ x_min  x_max], [yj yj], 'c');
    elseif x_max > 0
        plot(0, yj, 'c.');
    end
end
% Dibujar contornos del semianillo en negro
th = linspace(0, pi, 300);
plot(R_ext*cos(th), R_ext*sin(th), 'k', 'LineWidth', 1.5);
plot(R_int*cos(th), R_int*sin(th), 'k', 'LineWidth', 1.5);
plot([R_int R_ext]*cos(0),  [R_int R_ext]*sin(0),  'k', 'LineWidth', 1.5);
plot([R_int R_ext]*cos(pi), [R_int R_ext]*sin(pi), 'k', 'LineWidth', 1.5);
% Ajustar vista
xlim([-2.5 2.5]); ylim([-0.1 2.5]);
title('Mallado que representa los puntos interiores del sólido');
hold off;

2 Dibujar la temperatura del sólido

• La temperatura se incrementa a medida que ρ se aproxima a 1, porque el denominador del logaritmo disminuye.

• Esto indica un centro de calor concentrado en los puntos que se encuentran a una distancia de 1 desde el origen.

• La temperatura disminuye gradualmente al alejarnos del centro (hacia [math]\rho = 2[/math]).

Código MATLAB y representación:

Figura 2.1: Temperatura T(x,y).

clc; clear; close all;
% Malla estructural 2D con temperatura T(x,y) = (x - y)^2 en semianillo
% Parámetros de malla
r_min = 1;
r_max = 2;
paso = 0.1;
% Vectores de radios y ángulos
r = r_min:paso:r_max;
theta = linspace(0, pi, round(pi/paso)+1);  % asegura que incluya pi
% Preparar figura
figure;
hold on;
axis equal;
grid on;
title('Distribución de temperatura T(x,y) en la sección');
xlabel('X');
ylabel('Y');
% Pintar cada sector con color según temperatura
for i = 1:length(r)-1
   for j = 1:length(theta)-1
       r1 = r(i); r2 = r(i+1);
       th1 = theta(j); th2 = theta(j+1);
       % Coordenadas de los vértices del sector
       x = [r1*cos(th1), r2*cos(th1), r2*cos(th2), r1*cos(th2)];
       y = [r1*sin(th1), r2*sin(th1), r2*sin(th2), r1*sin(th2)];
       % Centro del sector para evaluar temperatura
       xc = mean(x);
       yc = mean(y);
       T = (xc - yc)^2;
       % Rellenar el sector con color proporcional a T
       fill(x, y, T, 'EdgeColor', 'k');   % borde negro
   end
end
% Dibujar líneas radiales
for j = 1:length(theta)
   plot(r.*cos(theta(j)), r.*sin(theta(j)), 'k-');
end
% Dibujar líneas circulares
for i = 1:length(r)
   plot(r(i)*cos(theta), r(i)*sin(theta), 'k-');
end
% Contornos más gruesos
th = linspace(0, pi, 300);
plot(r_max*cos(th), r_max*sin(th), 'k', 'LineWidth', 1.5); % arco exterior
plot(r_min*cos(th), r_min*sin(th), 'k', 'LineWidth', 1.5); % arco interior
plot([r_min r_max]*cos(0),  [r_min r_max]*sin(0),  'k', 'LineWidth', 1.5); % lateral derecho
plot([r_min r_max]*cos(pi), [r_min r_max]*sin(pi), 'k', 'LineWidth', 1.5); % lateral izquierdo
% Ajustes visuales
colormap(jet);     % mapa de colores tipo térmico
colorbar;          % barra de escala de temperatura

• Para ver la temperatura en cada punto del semianillo, se emplea un mapa de colores (colorbar).

• Las áreas de temperatura más alta (cerca de \rho = 1) están señaladas por los colores más fuertes.

• Las áreas más frías (alrededor de [math]\rho = 2[/math]) están representadas por los colores menos intensos.

3 Calcular el gradiente y dibujarlo como campo vectorial

Campo de desplazamiento en polares: [math]\vec{u}(\rho,\theta) = \frac{1}{5}(\rho - 1)\rho \,\vec{e}_{\rho}[/math]

Función de temperatura: [math]T(x,y) = (x - y)^2[/math]

Cálculo del gradiente:

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

[math]\frac{\partial T}{\partial x} = 2 - 2y;[/math] [math] \frac{\partial T}{\partial y} = -2 + 2x[/math]

[math]\nabla T(x,y) = (\,2 - 2y,\,-2 + 2x\,)[/math]

Curvas de nivel de T:

La condición de las curvas de nivel es:

[math]T(x,y) = c\;\;\Rightarrow\;\; (x - y)^2=c;[/math] [math](x - y)^2 = c\quad\Rightarrow\quad y=x\mp\sqrt{c}[/math]

* El dominio es ρ ∈ [1,2], θ ∈ [0,π], es decir, el semianillo del plano superior. Se debe a que la sección longitudinal del semianillo recorta las curvas de nivel y el campo.

* Las curvas de nivel son rectas que, al intersectar el semianillo generan segmentos rectos dentro de la figura.

Ortogonalidad de ∇T a las curvas de nivel:

Para ello, se considera el diferencial total de la función de temperatura:

[math]dT = \frac{\partial T}{\partial x}\,dx + \frac{\partial T}{\partial y}\,dy = \nabla T \cdot d\vec{r}[/math]

En una curva de nivel, donde [math]T(x,y) = c[/math], se cumple:

[math]dT = 0 \;\;\Rightarrow\;\; \nabla T \cdot d\vec{r} = 0[/math]

* Esto demuestra que el gradiente [math]\nabla T[/math] es ortogonal al vector tangente [math]d\vec{r}[/math] de la curva de nivel. Las flechas del gradiente son normales a cada curva de nivel porque, por construcción, el gradiente apunta en la dirección del máximo incremento de T y el valor de T es constante a lo largo de la curva.

Campo radial en la sección del arco y su conversión:

Campo dado en polares:

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

Campo en cartesianas:

[math]u_{\rho} = \frac{1}{5}(\rho - 1)\rho[/math]
[math]U_x = u_{\rho} \cos\theta;\quad U_y = u_{\rho} \sin\theta[/math]

* Dirección radial (sale del centro).

* Magnitud nula en [math]\rho = 1[/math] y máxima en [math]\rho = 2[/math] (dentro del arco mostrado).

Código MATLAB y representación

Figura 3.1: Campo del gradiente de la temperatura

clc; clear; close all;
% Paso deseado
paso = 0.1;
% Vectores con paso fijo
theta = 0:paso:pi;         % Ángulo de 0 a pi con paso 0.1
r = 1:paso:2;              % Radios de 1 a 2 con paso 0.1
% Malla polar
[Theta, R] = meshgrid(theta, r);
% Coordenadas cartesianas
X = R .* cos(Theta);
Y = R .* sin(Theta);
% Campo vectorial u(rho,theta) = (1/5)*(rho - 1)*rho * e_rho
U_rho = (1/5) * (R - 1) .* R;
U_x = U_rho .* cos(Theta);  % componente radial en x
U_y = U_rho .* sin(Theta);  % componente radial en y
% Función de temperatura T(x,y) = (x - y)^2
T = (X - Y).^2;
% Gradiente de T: ∇T = [∂T/∂x, ∂T/∂y]
dTdx = 2 - 2*Y;
dTdy = -2 + 2*X;
% Figura
figure;
hold on; axis equal; grid on;
% Dibujar malla del arco
for i = 1:length(theta)
    plot(X(:,i), Y(:,i), 'y');
end
for i = 1:length(r)
    plot(X(i,:), Y(i,:), 'y');
end
% Contornos del arco
th = linspace(0, pi, 300);
plot(2*cos(th), 2*sin(th), 'k', 'LineWidth', 1.5);  % Arco exterior
plot(1*cos(th), 1*sin(th), 'k', 'LineWidth', 1.5);  % Arco interior
plot([1 2]*cos(0),  [1 2]*sin(0),  'k', 'LineWidth', 1.5);  % Lateral 0
plot([1 2]*cos(pi), [1 2]*sin(pi), 'k', 'LineWidth', 1.5);  % Lateral pi
% Campo vectorial ∇T en rojo
quiver(X, Y, dTdx, dTdy, 'r', 'LineWidth', 1.2);
% Curvas de nivel de T en azul
contour(X, Y, T, 15, 'LineColor', 'b', 'LineWidth', 1.2);
% Campo u en verde
quiver(X, Y, U_x, U_y, 'g', 'LineWidth', 1.2);
% Ajustes finales
title('Campo ∇T (rojo), u (verde) y curvas de nivel T(x,y) (azul)');
xlabel('x'); ylabel('y');
xlim([-2.5 2.5]); ylim([-1.5 3]);
legend('∇T', 'Curvas de nivel T', 'Campo u');
hold off;

4 Dibujar el campo de vectores en los puntos del mallado del sólido.

Líneas coordenadas en cartesianas: [math]\vec{u}(x,y) = \frac{\cos(\pi y)}{10}\,\vec{i}[/math]

Código MATLAB y representación:

Figura 4.1: Campo del gradiente

% Medio arco: rho in [1,2], theta in [0,pi]
rho = linspace(1,2,20);
theta = linspace(0,pi,40);
[R,Theta] = meshgrid(rho,theta);
% Coordenadas cartesianas
X = R .* cos(Theta);
Y = R .* sin(Theta);
% Campo radial hacia fuera:

Esta sección ilustra la manera de dibujar un campo vectorial que irradia desde el centro en un área con forma de arco. Los vectores siempre señalan hacia afuera, y su tamaño se ajusta a cuán lejos están del punto central, lo cual demuestra una acción de crecimiento.

En la ilustración se aprecia sin problemas:

• Todos los vectores se dirigen hacia afuera, siguiendo la curva del arco.

• Las flechas se hacen más largas a medida que nos acercamos a la distancia radial de 2.

• El borde del arco parcial nos permite observar que el campo solo existe en esa zona.

5 Dibujar el sólido antes y después del desplazamiento

Código MATLAB y representación:

Figura 5.1:Semianillo en reposo

Figura 5.2:Semianillo desplazado por campo radial

clc; clear; close all;
% Parámetros del semianillo
R_ext = 2; % radio exterior
R_int = 1; % radio interior
n_ang = 40; % divisiones angulares
n_rad = 15; % divisiones radiales
% Ángulos y radios
theta = linspace(0, pi, n_ang);
r = linspace(R_int, R_ext, n_rad);
% Crear malla en coordenadas polares
[Theta, R] = meshgrid(theta, r);
% Coordenadas cartesianas en reposo
X0 = R .* cos(Theta);
Y0 = R .* sin(Theta);
% Campo de desplazamiento radial
U_r = (1/5) * (R - 1) .* R;
% Coordenadas desplazadas
R1 = R + U_r;
X1 = R1 .* cos(Theta);
Y1 = R1 .* sin(Theta);
% Dibujo con subplot
figure('Color','w');
% Subplot 1: semianillo en reposo
subplot(1,2,1);
hold on; axis equal; grid on;
title('Semianillo en reposo');
xlabel('x'); ylabel('y');
% Dibujar malla en amarillo
for i = 1:n_ang
plot(X0(:,i), Y0(:,i), 'b');
end
for i = 1:n_rad
plot(X0(i,:), Y0(i,:), 'b');
end
% Contornos en negro
th = linspace(0, pi, 300);
plot(R_ext*cos(th), R_ext*sin(th), 'k', 'LineWidth', 1.5); % arco exterior
plot(R_int*cos(th), R_int*sin(th), 'k', 'LineWidth', 1.5); % arco interior
plot([R_int R_ext]*cos(0), [R_int R_ext]*sin(0), 'k', 'LineWidth', 1.5); % lateral 1
plot([R_int R_ext]*cos(pi), [R_int R_ext]*sin(pi), 'k', 'LineWidth', 1.5); % lateral 2
xlim([-2.5 2.5]); ylim([-1.5 3]);
% Subplot 2: semianillo desplazado
subplot(1,2,2);
hold on; axis equal; grid on;
title('Semianillo desplazado por campo radial');
xlabel('x'); ylabel('y');
% Dibujar malla desplazada
for i = 1:n_ang
plot(X1(:,i), Y1(:,i), 'b');
end
for i = 1:n_rad
plot(X1(i,:), Y1(i,:), 'b');
end
% Contornos desplazados
th = linspace(0, pi, 300);
R_ext_d = R_ext + (1/5)*(R_ext - 1)*R_ext;
R_int_d = R_int + (1/5)*(R_int - 1)*R_int;
plot(R_ext_d*cos(th), R_ext_d*sin(th), 'k', 'LineWidth', 1.5); % arco exterior desplazado
plot(R_int_d*cos(th), R_int_d*sin(th), 'k', 'LineWidth', 1.5); % arco interior desplazado
% Líneas laterales desplazadas (cerrando por abajo)
x_lateral = [R_int_d, R_ext_d];
y_lateral = [0, 0]; % porque sin(0) = sin(pi) = 0
plot(x_lateral, y_lateral, 'k', 'LineWidth', 1.5); % lateral derecho (θ = 0)
plot(-x_lateral, y_lateral, 'k', 'LineWidth', 1.5); % lateral izquierdo (θ = π)
xlim([-2.5 2.5]); ylim([-1.5 3]);
hold off;

En esta sección se ilustra de manera gráfica cómo cambia el sólido cuando se le aplica el campo de desplazamiento.

Este campo señala que cada punto del sólido se mueve radialmente hacia afuera y que la distancia al centro determina la magnitud del desplazamiento: es cero en [math]\rho = 1[/math] y llega a su máximo en [math]\rho = 2[/math].

A fin de observar este fenómeno:

• El sólido original, que es un semianillo entre los radios 1 y 2, se dibuja primero.

• Después, se representa el sólido desplazado, en el que cada punto ha sido trasladado de acuerdo con el campo [math]\vec{u}[/math].

• Se pueden comparar de manera directa la geometría antes y después del desplazamiento, ya que ambos se encuentran en la misma figura usando subplot.

Este tipo de movimiento simula una expansión radial, como si el material se estuviera empujando desde el centro hacia afuera. Es característico en modelos de ondas sísmicas tipo S o en procedimientos de expansión térmica que se llevan a cabo en estructuras con forma circular.

6 Calcula la divergencia en todos los puntos del sólido

Código MATLAB y representación:

Figura 6:Divergencia

% Parámetros del dominio
rho_min = 1; rho_max = 2;
theta_min = 0; theta_max = pi;

% Resolución de la malla
Nr = 80;   % puntos en rho
Nt = 160;  % puntos en theta

% Malla polar
[theta, rho] = meshgrid(linspace(theta_min, theta_max, Nt), ...
                        linspace(rho_min, rho_max, Nr));

% Conversión a coordenadas cartesianas
x = rho .* cos(theta);
y = rho .* sin(theta);

% Divergencia: (1/5)(rho - 1)rho cos(theta)
div_u = (1/5) * (rho - 1) .* rho .* cos(theta);

% --- Gráfico ---
figure;

% Mapa de calor en azul
contourf(x, y, div_u, 50, 'LineColor', 'none');
colormap(parula);   % paleta azul-amarillo, puedes usar 'winter' para azul puro
colorbar;
title('Divergencia del campo (\nabla \cdot u)');
xlabel('x'); ylabel('y');
axis equal tight;

Se realiza el cálculo de la divergencia [math]\nabla \cdot \vec{u}[/math] en cada uno de los puntos del sólido y se representa gráficamente.

El campo de desplazamiento considerado es:

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

La divergencia es una medida del cambio de volumen local que ocurre a causa del desplazamiento. En este caso, nos permite identificar las zonas del semianillo donde el material se expande o se comprime.

Leyenda visual (colorbar en MATLAB):

La franja de colores que aparece en el gráfico es una escala visual que traduce los colores del mapa en valores numéricos de divergencia.

Interpretación gráfica:

• El valor de [math]\nabla \cdot \vec{u}[/math] se calcula en cada punto del semianillo.

• En el gráfico, ese valor se representa mediante un color.

• La barra lateral indica la relación entre el color y el número.

Por ejemplo:

🔵 Azul oscuro → valores negativos → el campo está "comprimido" localmente.

🔷 Azul claro → valores cercanos a cero → el volumen apenas cambia.

⚪ Colores más claros → valores positivos → el campo está "ampliando" el volumen local.

Interpretación física:

• Colores más intensos → mayor magnitud de divergencia → más fuerte el cambio en el volumen.

• Colores próximos a cero → poco cambio en el volumen.

De esta manera, puedes observar de un solo vistazo dónde el campo produce expansión (divergencia positiva) y dónde se genera compresión (divergencia negativa).

7 Calcular el rotacional en todos los puntos del sólido

Código MATLAB y representación:

Figura 4.2: Campo radial en medio arco.

% Medio arco: radios 1 a 2, ángulos 0 a pi
r1 = 1; r2 = 2; t1 = 0; t2 = pi;
dr = 0.05; dt = 0.05;
r = r1:dr:r2; t = t1:dt:t2;
[TT, RR] = meshgrid(t, r);
% Coordenadas cartesianas
X = RR .* cos(TT);
Y = RR .* sin(TT);
% Rotacional del campo radial: nulo en todo el dominio
C = zeros(size(RR));
% Puntos coloreados en azul
figure;
scatter(X(:), Y(:), 15, C(:), 'filled');
colormap('winter'); colorbar;   % 'winter' = gama azul
hold on;
% Contorno negro del medio arco
tt = linspace(t1, t2, 300);
plot(r1*cos(tt), r1*sin(tt), 'k', r2*cos(tt), r2*sin(tt), 'k');
plot([r1*cos(t1), r2*cos(t1)], [r1*sin(t1), r2*sin(t1)], 'k');
plot([r1*cos(t2), r2*cos(t2)], [r1*sin(t2), r2*sin(t2)], 'k');
axis equal;
xlabel('x'); ylabel('y');
title('campo radial en medio arco');
grid on; hold off;

Rotacional del campo radial en coordenadas polares

Considera el campo en el plano, en coordenadas polares:

[math]\vec{u}(\rho, \theta) = u_\rho(\rho, \theta) \, \vec{e}\rho + u\theta(\rho, \theta) \, \vec{e}_\theta[/math]

con

[math]u_\rho(\rho, \theta) = \frac{1}{\rho} (\rho - 1) \rho, \quad u_\theta(\rho, \theta) = 0[/math]

La componente [math]z[/math] del rotacional en 2D (plano) usando coordenadas polares es:

[math](\nabla \times \vec{u})z = \frac{1}{\rho} \frac{\partial}{\partial \rho} (\rho u\theta) - \frac{1}{\rho} \frac{\partial u_\rho}{\partial \theta}[/math]

Primer término:

• [math]\frac{1}{\rho} \frac{\partial}{\partial \rho} (\rho u_\theta) = \frac{1}{\rho} \frac{\partial}{\partial \rho} (\rho \cdot 0) = 0[/math]

Segundo término:

• [math]\frac{1}{\rho} \frac{\partial u_\rho}{\partial \theta} = \frac{1}{\rho} \frac{\partial}{\partial \theta} \left( \frac{1}{\rho} (\rho - 1) \rho \right) = 0[/math]

porque [math]u_\rho[/math] no depende de [math]\theta[/math].

Por tanto:

[math]\nabla \times \vec{u} = 0 \cdot \vec{k} = 0[/math]

¿Qué puntos sufren un mayor rotacional?

El cálculo del rotacional da cero en todo el dominio porque:

[math]u_\theta = 0[/math]

[math]u_\rho[/math] solo depende de [math]\rho[/math] y no de [math]\theta[/math]

Por tanto:

[math]\nabla \times \vec{u} = 0[/math] en todo el arco.

Ningún punto sufre rotación. El desplazamiento es solo hacia fuera del centro del arco. Por tanto, los puntos se expanden radialmente, pero no giran.

8 Tensor de deformaciones y tensiones normales en el medio elástico

Código MATLAB y representación:

Figura 8.1: Campo dirección radial

Figura 8.2: Campo dirección tangencial

clc; clear; close all;
% Parámetros geométricos
r_min = 1;
r_max = 2;
paso  = 0.1;
% Parámetros materiales
lambda = 1;
mu     = 1;
% Malla polar
r = r_min:paso:r_max;
theta = linspace(0, pi, round(pi/paso)+1);
[R, Theta] = meshgrid(r, theta);
% Coordenadas cartesianas
X = R .* cos(Theta);
Y = R .* sin(Theta);
% Componentes de deformación
a = (1/5) * (2*R - 1);   % derivada campo vectorial
b = (1/5) * (R - 1);    
div_u = (1/5) * (3*R - 2);  %divergencia
% Tensiones normales
sigma_rr = lambda .* div_u + 2 * mu .* a; %componente radial (efecto directo fuerza)
sigma_tt = lambda .* div_u + 2 * mu .* b; %componente tangencial (efecto indirecto fuerza)
% Vectores base en cada punto
e_r_x = cos(Theta);
e_r_y = sin(Theta);
e_t_x = -sin(Theta);
e_t_y = cos(Theta);
% Campo vectorial: flechas de tensión
U_rr_x = sigma_rr .* e_r_x;
U_rr_y = sigma_rr .* e_r_y;
U_tt_x = sigma_tt .* e_t_x ./ R;
U_tt_y = sigma_tt .* e_t_y ./ R;
% Visualización
figure('Position',[100 100 1000 500]);
% σ_rr como campo vectorial
subplot(1,2,1); hold on; axis equal;
grid on
title('\sigma_{\rho\rho}: campo en dirección radial');
xlabel('X'); ylabel('Y');
quiver(X, Y, U_rr_x, U_rr_y, 0.5, 'r');  % flechas rojas
plot_mallado(r, theta);
% --- σ_tt como campo vectorial ---
subplot(1,2,2); hold on; axis equal;
grid on
title('\sigma_{\theta\theta}: campo en dirección tangencial');
xlabel('X'); ylabel('Y');
quiver(X, Y, U_tt_x, U_tt_y, 0.5, 'b');  % flechas azules
plot_mallado(r, theta);
% Función auxiliar para mallado y contornos
function plot_mallado(r, theta)
   for j = 1:length(theta)
       plot(r.*cos(theta(j)), r.*sin(theta(j)), 'k-', 'LineWidth', 0.5);
   end
   for i = 1:length(r)
       plot(r(i)*cos(theta), r(i)*sin(theta), 'k-', 'LineWidth', 0.5);
   end
   th = linspace(0, pi, 300);
   plot(r(end)*cos(th), r(end)*sin(th), 'k', 'LineWidth', 1.5);
   plot(r(1)*cos(th), r(1)*sin(th), 'k', 'LineWidth', 1.5);
   plot([r(1) r(end)]*cos(0),  [r(1) r(end)]*sin(0),  'k', 'LineWidth', 1.5);
   plot([r(1) r(end)]*cos(pi), [r(1) r(end)]*sin(pi), 'k', 'LineWidth', 1.5);
end

Gradiente campo vectorial:

Como [math]\vec{u}[/math] es radial y solo depende de [math]\rho[/math], es un tensor simétrico y por tanto su parte simétrica es él mismo.

El gradiente del campo vectorial [math]\vec{u}(\rho, \theta)[/math] se expresa como:

[math]\nabla \vec{u} =\left(\begin{array}{cc}\frac{1}{5}(2\rho - 1) & 0 \\0 & \frac{1}{5}(\rho - 1)\end{array}\right)[/math]

El término inferior derecho se obtiene como:

[math]\frac{u_{\rho}}{\rho} = \frac{1}{\rho} \cdot \frac{1}{5}(\rho^2 - \rho) = \frac{1}{5}(\rho - 1)[/math]

Divergencia:

[math]\frac{\partial}{\partial \rho} (\rho\, u_{\rho}) = \frac{1}{5}(3\rho^2 - 2\rho)[/math]