Ver el código de «Mallado Arco 1 (grupo 59)»

{{ TrabajoED | Arco1. Grupo 59 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC25/26|2025-26]] | 
*Guillermo Pineros Quero
*Ignacio Sueiras Oviedo
*Francisco José Saiz Cebrián
*Álvaro Villar Rosado
*Javier Ruiz Sáenz de Jubera}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC25/26]]
==Introducción==
"En este artículo se examina el proceso de parametrización y mallado de una placa plana 2D con forma de arco circular (anillo), delimitada por un radio interior de 1 y uno exterior de 2. Definimos dos variables físicas sobre esta geometría: la temperatura <math>T(x,y)</math> en coordenadas cartesianas y un vector de desplazamientos <math>\vec{u}(\rho,\theta)</math> en coordenadas cilíndricas. La temperatura se rige por la función      <math>T(x,y)=(x-y)^2</math> , mientras que los desplazamientos obedecen a <math>\vec{u}</math> = <math>\frac{1}{5} (\rho-1)\rho^2 \sin\theta(\vec{e}_\theta)</math>. El objetivo principal es formalizar la descripción del dominio y sus funciones para facilitar el análisis numérico mediante discretización. El texto se complementa con códigos de MATLAB, fotografías y representaciones gráficas en cada apartado."
=== Póster del trabajo ===
A continuación subimos el enlace de descarga del póster del trabajo.

-[[Archivo:Arco_1_Grupo59.pdf|thumb|left|Póster del artículo resumido]]

== Mallado del arco ==
El mallado de un arco consiste en dividir su superficie en puntos organizados, lo que permite describir su geometría y analizar cómo cambian diferentes magnitudes dentro de él.
El arco está comprendido entre los radios de valor 1 y 2, que al pasarlo a coordenadas polares la región estará definida por rho(ρ) y theta(θ), y sus valores [1,2]x[0,π]

=== Representación del mallado ===
[[Archivo:mallado.del.arco.png|500px|thumb|right|]]
{{matlab|codigo=
% Mallado del arco alineado

% Elegimos divisiones que garanticen que:
% θ = 0, pi/2, pi están incluidos 
% r incluye exactamente 1 y 2    

theta = linspace(0, pi, 20);   % mallado moderado, incluye 0, pi/2 y pi
r = linspace(1, 2, 15);        % radios exactos 1 y 2

[R, TH] = meshgrid(r, theta);

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

figure;
hold on;

col = [0.3 0.6 1];  % azul clarito

% Mallado tipo red (más grueso) 
plot(X,  Y, 'Color', col, 'LineWidth', 1.1);   % líneas en theta
plot(X', Y', 'Color', col, 'LineWidth', 1.1);  % líneas en r

% Contornos superior interior y exterior 
theta_cont = linspace(0, pi, 400);

plot(1*cos(theta_cont), 1*sin(theta_cont), 'Color', col, 'LineWidth', 2);
plot(2*cos(theta_cont), 2*sin(theta_cont), 'Color', col, 'LineWidth', 2);

% Bases 
plot([-1 -2], [0 0], 'Color', col, 'LineWidth', 2);  % base izquierda
plot([ 1  2], [0 0], 'Color', col, 'LineWidth', 2);  % base derecha


% Ajustes
axis equal
grid on
title('Mallado del arco')
xlabel('x')
ylabel('y')
hold off;
}}

==  Temperatura en el arco  ==
La temperatura del sólido proviene de un foco de calor muy concentrado en los puntos
que están a distancia 1 del origen. La función es: '''<big>𝑇(𝑥, 𝑦) = (𝑥 − 𝑦)<sup>2</sup></big>'''

=== Código ===
Este código crea una imagen que representa como afecta la temperatura al arco, en donde se ven los distintos valores de temperatura representados mediante colores, donde el azul representa las temperaturas mas bajas y el rojo representa las temperaturas mas altas.
[[Archivo:funciontemperatura.png|500px|thumb|right|Representación de la temperatura en el arco]]
{{matlab|codigo=
% Mallado del arco
theta = linspace(0, pi, 200);   % más denso para curva suave
r = linspace(1, 2, 200);        % radio interior 1, exterior 2
[R, TH] = meshgrid(r, theta);

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

% Funcion temperatura
T = (X - Y).^2;

% Representación 2D
figure;
pcolor(X, Y, T);     % representación 2D en colores
shading interp       % suavizado del color
colormap(jet)
colorbar;

title('Temperatura T(x,y) = (x - y)^2 sobre el arco')
xlabel('x')
ylabel('y')
axis equal           % mantiene proporciones reales
hold on

%% Contorno negro del arco 
theta_cont = linspace(0, pi, 400);

% Borde interior
plot(1*cos(theta_cont), 1*sin(theta_cont), 'k-', 'LineWidth', 1);

% Borde exterior
plot(2*cos(theta_cont), 2*sin(theta_cont), 'k-', 'LineWidth', 1);

% Radio izquierdo
plot([1 2]*cos(0), [1 2]*sin(0), 'k-', 'LineWidth', 1);

% Radio derecho
plot([1 2]*cos(pi), [1 2]*sin(pi), 'k-', 'LineWidth', 1);

hold off
}}

== Gradiente de temperatura ==
===Cálculo del gradiente===
El gradiente de una función escalar (∇f) es un campo vectorial cuyas componentes son las derivadas parciales de la función. Geométricamente, el vector gradiente apunta en la dirección de máximo crecimiento de la función y es ortogonal a sus curvas de nivel.

La fórmula para calcularlo es: 

∇f(x,y) = ( ∂f/∂x , ∂f/∂y )

Como la temperatura viene dada por la siguiente fórmula: 

T(x, y) = (x - y)²

El gradiente de temperatura será:

∇T = (2(x - y), -2(x - y))
===Representación gráfica del gradiente===
[[Archivo:Gradiente apartado 3.png|500px|thumb|right|Gradiente de temperatura]]
{{matlab|codigo=
%% Representación del gradiente
clc; clear; close all;

% Dominio 
rho_min = 1;  rho_max = 2;
th_min = 0;   th_max = pi;

[Rho, Theta] = meshgrid(linspace(rho_min, rho_max, 100), ...
                        linspace(th_min, th_max, 200));

% Convertimos a Cartesianas
X = Rho .* cos(Theta);
Y = Rho .* sin(Theta);

%% 2. Definición de Temperatura y Gradiente
% Temperatura T = (x - y)^2
T_val = (X - Y).^2;

% Cálculo del gradiente
Grad_X = 2 .* (X - Y);
Grad_Y = -2 .* (X - Y);

% Representación gráfica
figure('Name', 'Trabajo M - Pregunta 3 (Estilo Suave)', 'Color', 'w');
hold on;

contourf(X, Y, T_val, 50, 'LineStyle', 'none');

% Creamos la gráfica
N = 256;
r = linspace(0.2, 1, N)'; 
g = linspace(0.5, 1, N)';
b = linspace(1, 0, N)';   
mi_mapa = [r, g, b];      % Matriz de colores

colormap(mi_mapa);
c = colorbar;
c.Label.String = 'Temperatura T = (x-y)^2';

% Cantida de flechas
step_r = 8;  % Salta cada 8 puntos en el radio
step_t = 12; % Salta cada 12 puntos en el ángulo

idx_r = 1:step_r:size(X,1);
idx_t = 1:step_t:size(X,2);

% Dibujo flechas
quiver(X(idx_r, idx_t), Y(idx_r, idx_t), ...
       Grad_X(idx_r, idx_t), Grad_Y(idx_r, idx_t), ...
       'k', ...                 % Color negro flechas
       'LineWidth', 1.0, ...    % Grosor flechas
       'AutoScaleFactor', 0.8); % Tamaño flechas

axis equal;
xlim([-2.5 2.5]); ylim([0 2.5]); % Márgenes 
title('Gradiente sobre T');
xlabel('x'); ylabel('y');

box on; 
hold off;
}}
== Campo de Vectores ==
El campo vectorial viene dado por la siguiente fórmula: 
<div style="text-align:center;">
<math>\vec{u}(\rho,\theta)=\frac{1}{5}(\rho - 1)\rho^{2}\sin\theta\,\vec{e}_{\theta}</math>
</div>

=== código ===
[[Archivo:apartado.4.png|500px|thumb|right|Campo de vectores]]
{{matlab|codigo=
clear; clc; close all;

%% --- 1. Definir el mallado polar del sólido ---
h = 0.1;
rho_min = 1;
rho_max = 2;
theta_min = 0;
theta_max = pi;

rho_vec   = rho_min:h:rho_max;
theta_vec = theta_min:h:theta_max;

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

% Convertir la malla a cartesianas (puntos del sólido)
X = R .* cos(Th);
Y = R .* sin(Th);

%% --- 2. Definir el campo vectorial u en e_theta ---
U_theta = (1/5) * (R - 1) .* R.^2 .* sin(Th);

% Componentes de los vectores unitarios e_theta en cartesianas
e_th_x = -sin(Th);
e_th_y =  cos(Th);

% Convertir el campo a componentes cartesianas
U_x = U_theta .* e_th_x;
U_y = U_theta .* e_th_y;

%% --- 3. Dibujar SOLO el campo vectorial u en los puntos de la malla ---
figure;
hold on;
axis equal;
grid off;
set(gcf,'Color','w');

% Dibujar las flechas del campo vectorial u
quiver(X, Y, U_x, U_y, 'b', 'LineWidth', 1.2, 'MaxHeadSize', 0.6);

%% --- 4. Dibujar el borde del sólido con UNA línea negra ---
% Lado 1: arco exterior
plot(rho_max*cos(theta_vec), rho_max*sin(theta_vec), 'k', 'LineWidth', 2);

% Lado 2: arco interior
plot(rho_min*cos(theta_vec), rho_min*sin(theta_vec), 'k', 'LineWidth', 2);

% Lado 3: cierre izquierdo (theta = pi)
plot([rho_min*cos(pi), rho_max*cos(pi)], [rho_min*sin(pi), rho_max*sin(pi)], 'k', 'LineWidth', 2);

% Lado 4: cierre derecho (theta = 0)
plot([rho_min*cos(0), rho_max*cos(0)], [rho_min*sin(0), rho_max*sin(0)], 'k', 'LineWidth', 2);

title('Campo \vec u sobre la malla del sólido (bordeado en negro)');
xlabel('x'); ylabel('y');
hold off;
}}

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

===Código===
[[Archivo:grafica_1.png|500px|thumb|right|Gráfica 1]]
[[Archivo:Grafica_2.png|500px|thumb|right|Gráfica 2]]
[[Archivo:apartado.5.png|500px|thumb|right|Gráfica 3]]
{{matlab|codigo=
%% 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
}}

==Divergencia==
[[Archivo:apartado.6.png|500px|thumb|right|Divergencia]]
===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>

===Código y representación===
{{matlab|codigo=
%% 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(redbluecmap); 

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

}}

== Rotacional ==
El rotacional de un campo vectorial <math>|\nabla \times \vec{u}|</math> mide cuánto y en qué dirección tiende a girar el campo, es decir, si el campo hace que una partícula colocada en él empiece a rotar. Si el rotacional es distinto de 0 entonces la partícula rotará, en cambio si el rotacional es igual a 0 entonces la partícula no rotará. 
La fórmula del campo gradiente es <center><math>\nabla×\vec u(ρ,θ) = \frac{1}{\rho}\begin{vmatrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ  & \vec v_z \end{vmatrix}</math></center> 
El rotacional de un campo vectorial mide su tendencia a rotar en un punto.

La formula desarrollada queda de la siguiente forma:
<center><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></center>
<br />
=== Código ===
[[Archivo:apartado.7.png|500px|thumb|right|Rotacional]]
{{matlab|codigo=
clear; clc; close all;

% Parámetros del sólido (anillo semicircular)
h = 0.02;
rho_min = 1;
rho_max = 2;
theta_min = 0;
theta_max = pi;

rho_vec   = rho_min:h:rho_max;
theta_vec = theta_min:h:theta_max;
[R, Th] = meshgrid(rho_vec, theta_vec);

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

% Magnitud del rotacional
curl_z = (1/5) .* R .* (4.*R - 3) .* sin(Th);
curl_abs = abs(curl_z);

% --- Visualización ---
figure;
hold on; axis equal; grid off;
set(gcf,'Color','w');

contourf(X, Y, curl_abs, 30, 'LineStyle','none');
colorbar;
colormap(parula);

title('|∇×u| en el sólido');
xlabel('x'); ylabel('y');

% --- Marcar máximo explícitamente en (0,2) ---
xmax = 0;
ymax = 2;
plot(xmax, ymax, 'r*', 'MarkerSize', 14, 'LineWidth', 2);
text(xmax+0.05, ymax, '  Máximo rotacional (0,2)', 'Color','r','FontSize',10);

hold off;
}}
== Tensiones normales respecto a los planos tangencial y radial ==
=== Explicación de los cálculos ===
El cálculo de las tensiones normales en el arco se basa en la Ley de Hooke para un material elástico lineal que se comporta igual en todas las direcciones, 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{\epsilon} = \frac{1}{2}(\nabla \vec{u} + (\nabla \vec{u})^T)</math>, el tensor de deformaciones, que indica cómo cambia localmente la forma del material.

<math>\mathbf{I}</math>, la matriz identidad, que aplica uniformemente el efecto del cambio volumétrico.

Los coeficientes de Lamé:

<math>\lambda</math>, que mide la resistencia a cambios de volumen.

<math>\mu</math>, que mide la resistencia al corte (cizalladura).

Para simplificar, se toma <math>\lambda = 1</math> y <math>\mu = 1</math>, de modo que la Ley de Hooke se reduce a:
<math>\mathbf{\sigma = (\nabla \cdot \vec{u}) + 2\epsilon}</math>

Campo de desplazamientos en coordenadas polares:
<math>\vec{u}(\rho,\theta) = 0,\vec{e}{\rho} + \frac{1}{5}(\rho-1)\rho^2 \sin\theta\vec{e}{\theta}</math>
Es decir, el desplazamiento es solo tangencial y no hay componente radial.

Componentes del tensor de deformaciones (enumeradas):

<math>\epsilon_{\rho\rho} = \frac{\partial u_\rho}{\partial \rho} = 0</math>

<math>\epsilon_{\theta\theta} = \frac{1}{\rho} \frac{\partial u_\theta}{\partial \theta} + \frac{u_\rho}{\rho} = \frac{1}{5} (\rho^2 - \rho) \cos\theta</math>

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

La traza del tensor de deformaciones, <math>\mathrm{tr}(\epsilon) = \epsilon_{\rho\rho} + \epsilon_{\theta\theta}</math>, indica el cambio volumétrico relativo en el punto del material:
<math>\mathrm{tr}(\epsilon) = 0 + \frac{1}{5} (\rho^2 - \rho) \cos\theta = \frac{1}{5} (\rho^2 - \rho) \cos\theta</math>

Tensiones normales usando la Ley de Hooke simplificada (
𝜆
=
𝜇
=
1
λ=μ=1):

<div style="text-align:center;"> <math>\mathbf{\sigma_{\rho\rho}} = \vec{e}_{\rho} \cdot \mathbf{\sigma} \cdot \vec{e}_{\rho} = (\nabla \cdot \vec{u}) + 2\epsilon_{\rho\rho}</math> </div> <div style="text-align:center;"> <math>\mathbf{\sigma_{\rho\rho}} = \frac{1}{5} (\rho^2 - \rho) \cos\theta</math> </div> <div style="text-align:center;"> <math>\mathbf{\sigma_{\theta\theta}} = \vec{e}_{\theta} \cdot \mathbf{\sigma} \cdot \vec{e}_{\theta} = (\nabla \cdot \vec{u}) + 2\epsilon_{\theta\theta}</math> </div> <div style="text-align:center;"> <math>\mathbf{\sigma_{\theta\theta}} = \frac{3}{5} (\rho^2 - \rho) \cos\theta</math> </div>
=== Código y representaciones ===
[[Archivo:untitled8_1.png|500px|thumb|right|]]
[[Archivo:untitled8_2.png|500px|thumb|right|]]
{{matlab|codigo=
% Apartado 8 - tensiones normales sobre arco con ambas bases

% Mallado
theta = linspace(0, pi, 300);
rho   = linspace(1, 2, 200);
[TH, RHO] = meshgrid(theta, rho);

% Convertir a cartesianas para graficar
X = RHO .* cos(TH);
Y = RHO .* sin(TH);

% Campo dado: u_rho = 0, u_theta = (1/5)*(rho-1)*rho^2*sin(theta)
U_rho = zeros(size(RHO));
U_theta = (1/5) .* (RHO - 1) .* (RHO.^2) .* sin(TH);

% Pasos
dtheta = theta(2) - theta(1);
drho   = rho(2)   - rho(1);

% Componentes del tensor simétrico de deformaciones
[Utheta_dtheta, Utheta_drho] = gradient(U_theta, dtheta, drho);
[Urho_dtheta,   Urho_drho]   = gradient(U_rho,   dtheta, drho);

eps_rr = Urho_drho;
eps_tt = (1./RHO).*Utheta_dtheta + U_rho ./ RHO;
eps_rt = 0.5 .* ( (1./RHO).*Urho_dtheta + Utheta_drho - U_theta ./ RHO );

% Traza
trace_eps = eps_rr + eps_tt;

% Coeficientes de Lamé
lambda = 1; mu = 1;

% Tensiones
sigma_rr = lambda .* trace_eps + 2*mu .* eps_rr;
sigma_tt = lambda .* trace_eps + 2*mu .* eps_tt;

% Sustitución pedida: j = (1/rho)e_theta
sigma_jj_scaled = (1 ./ (RHO.^2)) .* sigma_tt;

% Línea base derecha (theta = 0) y base izquierda (theta = pi)
base_rho = linspace(1,2,200);
base_x_right = base_rho .* cos(0);    base_y_right = base_rho .* sin(0);
base_x_left  = base_rho .* cos(pi);   base_y_left  = base_rho .* sin(pi);

tc = linspace(0,pi,400);

%% Dibujo de tensiones normales radiales (2D)
figure;
pcolor(X, Y, sigma_rr);
shading interp; colormap(jet);
cb = colorbar;
cb.Label.String = 'Tensión (Pa)';
title('\sigma_{\rho\rho}  (tensión normal radial)');
xlabel('x'); ylabel('y'); axis equal; hold on;

% Contorno del arco
plot(1*cos(tc), 1*sin(tc), 'k-', 'LineWidth', 0.6);
plot(2*cos(tc), 2*sin(tc), 'k-', 'LineWidth', 0.6);
plot(base_x_right, base_y_right, 'k-', 'LineWidth', 0.6);   % base derecha (theta=0)
plot(base_x_left,  base_y_left,  'k-', 'LineWidth', 0.6);   % base izquierda (theta=pi)

hold off;

%% Dibujo de tensiones normales tangenciales (2D)
figure;
pcolor(X, Y, sigma_jj_scaled);
shading interp; colormap(jet);
cb = colorbar;
cb.Label.String = 'Tensión (Pa)';
title('\sigma_{\theta\theta}  (tensión normal tangencial)');
xlabel('x'); ylabel('y'); axis equal; hold on;

% Contorno del arco
plot(1*cos(tc), 1*sin(tc), 'k-', 'LineWidth', 0.6);
plot(2*cos(tc), 2*sin(tc), 'k-', 'LineWidth', 0.6);
plot(base_x_right, base_y_right, 'k-', 'LineWidth', 0.6);   % base derecha (theta=0)
plot(base_x_left,  base_y_left,  'k-', 'LineWidth', 0.6);   % base izquierda (theta=pi)

hold off;
}}

== Tensiones tangenciales respecto al plano ortogonal <math>\vec{e}_{\rho}</math> ==

Se calculará el punto donde el material sufre mayor cizalladura. 

La fórmula es la siguiente: <math>\mathbf{\tau} = \mathbf{\sigma} \cdot \vec{e}_{\rho} - (\vec{e}_{\rho} \cdot \mathbf{\sigma} \cdot \vec{e}_{\rho})\vec{e}_{\rho}</math>

La fórmula que debes usar para el vector de tensión tangencial (cizalladura) es: <math>|\mathbf{\tau}| = |\mathbf{\sigma_{\rho\theta}}|</math>

Y al aplicar la Ley de Hooke <math>\mathbf{\sigma_{\rho\theta} = 2 \epsilon_{\rho\theta}}</math> y los componentes de deformación, el resultado es:
=== Código ===

<math>|\mathbf{\tau}| = \left| \frac{1}{5} (2\rho^2 - \rho) \sin\theta \right|</math>
[[Archivo:apartado_9.jpg|500px|thumb|right|Tensiones tangenciales]]
{{matlab|codigo=
%% TENSIONES TANGENCIALES
clear; clc; close all;

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

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

% --- 2. CÁLCULO DE LA TENSIÓN ---
Tau = (1/5) * (2*R.^2 - R) .* sin(Th);
Tau_Mag = abs(Tau);

% --- 3. VISUALIZACIÓN ---
figure(9); clf; hold on; axis equal;
set(gcf, 'Color', 'w');
title('Tensión tangencial');
xlabel('x'); ylabel('y');

contourf(X, Y, Tau_Mag, 50, 'LineStyle', 'none');
c = colorbar;
ylabel(c,'Esfuerzo de Corte (Pa)');

% --- 4. Colormap amplio AZUL → ROJO (256 tonos) ---
n = 256;
azul = [0 0 1];       % azul puro
rojo = [1 0 0];       % rojo puro
cmap = [linspace(azul(1), rojo(1), n)', ...
        linspace(azul(2), rojo(2), n)', ...
        linspace(azul(3), rojo(3), n)'];
colormap(cmap);

% --- 5. Bordes exteriores en negro ---
t_border = 0:0.001:pi;
plot(2*cos(t_border),2*sin(t_border),'k','LineWidth',2);
plot(1*cos(t_border),1*sin(t_border),'k','LineWidth',2);
plot([-2,-1],[0,0],'k','LineWidth',2);
plot([1,2],[0,0],'k','LineWidth',2);

% --- 6. Punto de máxima tensión ---
max_val = max(Tau_Mag(:));
[fr, fc] = find(Tau_Mag == max_val, 1);
plot(X(fr,fc),Y(fr,fc),'wx','LineWidth',2,'MarkerSize',10);
text(X(fr,fc),Y(fr,fc)+0.1,'Máx','Color','w','FontWeight','bold');

axis([-2.2 2.2 0 2.2]);
grid off;
hold off;

}}

== 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| \frac{1}{\rho} \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| \frac{1}{\rho} \cdot \left[ \frac{1}{5} (2\rho^2 - \rho) \sin\theta \right] \right|</math>

[[Archivo:apartado.10.png|500px|thumb|right|Tensiones tangenciales]]
{{matlab|codigo=
clear; clc; close all;

%% --- 1. Mallado polar del sólido (corona semicircular) ---
h = 0.02;                 % paso fino
rho_min = 1;
rho_max = 2;
theta_min = 0;
theta_max = pi;

rho_vec = rho_min:h:rho_max;
theta_vec = theta_min:h:theta_max;

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

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

%% --- 2. Campo u = u_theta * e_theta ---
u_theta = (1/5) * (R - 1) .* R.^2 .* sin(Th);

% Cálculo del rotacional (componente z en 2D, analítico)
curl_z = (1/5) .* R .* (4.*R - 3) .* sin(Th);
curl_abs = abs(curl_z);

%% --- 3. Gráfica del rotacional (sin flechas) ---
figure;
hold on;
axis equal;
grid off;
set(gcf,'Color','w');

% Mapa de calor del rotacional
contourf(X, Y, curl_abs, 40, 'LineStyle','none');
colorbar;
clim([0 2]);      % forzamos rango para que el máximo destaque bien
colormap turbo;

title('Tension tangencial apartado 10');
xlabel('x'); ylabel('y');

% Borde del sólido
plot(rho_max*cos(theta_vec), rho_max*sin(theta_vec), 'k','LineWidth',1.8);
plot(rho_min*cos(theta_vec), rho_min*sin(theta_vec), 'k','LineWidth',1.8);
plot([-rho_max -rho_min], [0 0], 'k','LineWidth',1.8);
plot([rho_min rho_max], [0 0], 'k','LineWidth',1.8);

%  Marcar punto de máximo explícitamente en (0,2) en blanco
xmax = 0;
ymax = 2;
plot(xmax, ymax, 'w*', 'MarkerSize', 15, 'LineWidth', 2);
text(xmax+0.05, ymax-0.1, '  Máximo', 'Color','w','FontSize',11,'FontWeight','bold');

hold off;

%% --- 4. Mostrar valor numérico máximo en consola ---
fprintf('Máximo de |∇×u| = %.5g en (0,2)\n', curl_abs(find(X==0 & abs(Y-2)<1e-6, 1)));

disp('El punto con mayor rotacional es (0,2)');
disp('Porque |∇×u| crece con ρ y es máximo cuando sinθ = 1 (θ = π/2) y ρ = 2');
}}

==Cálculo de la masa aproximando la integral numéricamente==
La masa se calcula con la siguiente integral: M = ∫∫ d(ρ, θ) dA


El dominio es [1,2]x[0,pi] y por tanto los límites de integración serán:

Límete de rho: 1 ≤ ρ ≤ 2

Límite de theta: 0 ≤ θ ≤ π


El diferencial de área (Jacobiano) es: dA = ρ dρ dθ

La densidad viene dada por: d(ρ, θ) = 1 + e^(ρ² cos θ)


Por último la masa se calculará de la siguiente manera: 

M = ∫(de 0 a π) ∫(de 1 a 2) [1 + e^(ρ² cos θ)] · ρ dρ dθ = 26.26 u²

== Interpretación del trabajo ==

En el Trabajo K se modela una onda transversal en una placa rectangular mediante el campo de desplazamientos:

<math>\vec{u}(x,y,t)=\vec{a}\cos(\vec{b}\cdot\vec{r}_0 - c t)</math> con:

<math>\vec{a}=\frac{1}{10}\vec{i}</math>
<math>\vec{b}=\pi\vec{j}</math>
<math>t=0</math>

Esto se puede interpretar como una onda de corte (onda S) que se propaga en un medio elástico, como la corteza terrestre durante un terremoto.
Las ondas S (secundarias o de corte) son ondas transversales, en las que las partículas del material vibran perpendicularmente a la dirección de propagación de la onda.

En este caso, la dirección de propagación es la del vector:

<math>\vec{b}=\pi\vec{j}</math>

(eje 𝑦 y) y la vibración es en la dirección:<math>\vec{i}</math>(eje 𝑥 x), lo cual es típico de ondas S con polarización horizontal (onda SH).
Aplicación al Trabajo M (arco entre radios 1 y 2)
En el Trabajo M, el dominio es un sector de arco en coordenadas polares(𝜌,𝜃)(ρ,θ), con:<math>\rho\in[1,2]</math> y el campo de desplazamientos es:

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

Este campo representa un desplazamiento puramente tangencial: <math>\vec{e}_\theta</math> 

que también corresponde a una onda de corte pero en geometría curvilínea.
En sismología, las ondas S pueden viajar a lo largo de estructuras curvas en la corteza, como capas geológicas curvadas o interfaces entre materiales.

El factor:<math>(\rho-1)\rho^{2}\sin\theta</math> indica que el desplazamiento es máximo en la parte exterior del arco: <math>\rho = 2</math> y varía sinusoidalmente con:<math>\theta</math> lo que podría modelar una onda estacionaria o forzada en una placa curvada.

Relación con la temperatura

En ambos trabajos se incluye un campo de temperatura: <math>T(x,y)=(x-y)^2</math> que en el Trabajo M debe expresarse en coordenadas polares:<math>x=\rho\cos\theta</math>
<math>y=\rho\sin\theta</math>

Sustituyendo: <math>T(\rho,\theta)=\rho^{2}(\cos\theta - \sin\theta)^2</math>

Esta temperatura podría representar gradientes térmicos en la corteza debidos a fuentes geotérmicas, que afectan a las propiedades elásticas del medio y, por tanto, a la propagación de ondas sísmicas.

Conclusión

El Trabajo M puede interpretarse como el estudio de ondas de corte (S) en una región curvilínea de la corteza terrestre, como un arco de falla o una capa geológica curvada, donde los desplazamientos son transversales a la dirección radial.

La inclusión del campo de temperatura:<math>T(\rho,\theta)=\rho^{2}(\cos\theta - \sin\theta)^2</math> permite analizar cómo los gradientes térmicos influyen en la deformación y propagación de ondas en estructuras geológicas reales, algo relevante en sismología y geofísica.