Placa Plana (Grupo 20)

Trabajo realizado por estudiantes
Título	Placa Plana. Grupo 20
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Alejandro Trejo, Marcos Rodríguez-Barbero, Gloria García , Manuel Riesgo, Ángel De Miguel Pérez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Una placa plana es una superficie en dos dimensiones, que en este caso representa la sección longitudinal de una viga voladiza rectangular que está en una situación de voladizo apoyada sobre la pared izquierda, y ocupa la región [math](x, y) \in [0, 4] \times [f(g), g(x)][/math].
Las funciones [math]f(x) g(x)[/math] vienen definidas como [math]f(x) = \frac{x}{8}, \quad g(x) = 2 - \frac{x}{8}[/math].

Esta placa esta sometida a una temperatura que viene definida por la función [math]T(x, y) = (1 + (y - 1)^2)(4 - x)[/math].
Además, los puntos de la placa sufren un desplazamiento causado por una fuerza de magnitud desconocida. Este desplazamiento es tal, que partiendo de un vector posición genérico para cualquier punto de la placa ([math]\vec{r}_0(x, y) = x\vec{i} + y\vec{j}[/math]), al final del mismo el vector que nos sirve para determinar la posición de cualquier punto en la gráfica cambia a este:[math]\vec{r}_d(x, y) = \vec{r}_0(x, y) + \vec{u}(x, y)[/math].

Para realizar los diferentes cálculos, vamos a suponer que esta fuerza ha desplazado los puntos siguiendo el vector: [math]\vec{u}(\rho, \theta) = -\frac{1}{20}\rho^2 \cos \theta \vec{e}_\theta[/math].

2 Mallado de placa

Para empezar, vamos a realizar un mallado de la sección longitudinal que represente los puntos interiores del sólido. Hemos utilizado un paso de muestreo h=1/10, y establecido los ejes en el rectángulo acotado por: [math](x, y) \in [-1; 5] \times [-1; 3][/math].

% 1. Configuración
h = 1/10;
color_malla = [0, 0.6, 0.6]; % Color Cían Oscuro

% 2. Generar Coordenadas Paramétricas
% u: avanza a lo largo de la placa (0 a 4)
% v: avanza de abajo a arriba (0 a 2)
u = 0 : h : 4;
v = 0 : h : 2;
[U, V] = meshgrid(u, v);

% 3. Transformación Matemática
% Calculamos los bordes superior e inferior para cada punto X
y_abajo = U ./ 8;
y_arriba = 2 - (U ./ 8);

% Interpolamos: calculamos la Y real basada en la altura V
factor_altura = V ./ 2; % Va de 0 (abajo) a 1 (arriba)

X = U;
Y = y_abajo + factor_altura .* (y_arriba - y_abajo);

% 4. Visualización
figure('Color', 'w'); hold on;

% A) Dibujar las líneas VERTICALES
plot(X, Y, 'Color', color_malla, 'LineWidth', 0.5);

% B) Dibujar las líneas HORIZONTALES
plot(X', Y', 'Color', color_malla, 'LineWidth', 0.5);

% C) Dibujar el contorno 
plot(u, u./8, 'k-', 'LineWidth', 2); % Borde Abajo
plot(u, 2 - u./8, 'k-', 'LineWidth', 2); % Borde Arriba
plot([0 0], [0 2], 'k-', 'LineWidth', 2); % Borde Izquierdo
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); % Borde Derecho

% 5. Configuración final
axis([-1 5 -1 3]); % Zoom/Encuadre exacto
xlabel('x'); ylabel('y');
title('Mallado de la placa plana');
grid on;
box on;
hold off;

3 Curvas de nivel

% 1. Definición del Dominio y Mallado
h = 0.2; % Paso de muestreo
x = 0 : h : 4;
y = 0 : h : 2;
[X, Y] = meshgrid(x, y);

% 2. Definir la Geometría (Tu placa plana trapezoidal)
y_abajo = X ./ 8;
y_arriba = 2 - (X ./ 8);

% Creamos la máscara lógica (1 dentro, 0 fuera)
dentro = (Y >= y_abajo) & (Y <= y_arriba);

% 3. Definir la TEMPERATURA T(x,y)
% FÓRMULA ACTUALIZADA: T(x,y) = (1 + (y-1)^2)(4-x)
% Usamos .^ y .* para operaciones elemento a elemento
T_func = @(x,y) (1 + (y - 1).^2) .* (4 - x);

% Calculamos la temperatura en toda la malla
T = T_func(X, Y);

% 4. Aplicar la Máscara (Recortar la figura)
% Asignamos NaN a los puntos de fuera para que no se pinten
T(~dentro) = NaN;
% (Opcional) Recortar X e Y también ayuda a limpiar bordes en contour
X_plot = X; Y_plot = Y;
X_plot(~dentro) = NaN;
Y_plot(~dentro) = NaN;

% 5. Cálculo del Gradiente (Nabla T)
% gradient calcula las derivadas parciales numéricas
[dTdx, dTdy] = gradient(T, h, h);

% 6. Visualización
figure('Color', 'w', 'Position', [100 100 800 600]);
hold on;

% A) Pintar Curvas de Nivel (Isotermas)
% Usamos contour para ver las líneas de temperatura constante
[C, h_cont] = contour(X, Y, T, 20, 'LineWidth', 2);
clabel(C, h_cont, 'FontSize', 9, 'Color', 'k', 'FontWeight', 'bold');
colormap('jet');

% B) Pintar Campo Vectorial (Gradiente)
% quiver dibuja las flechas sobre la superficie
quiver(X, Y, dTdx, dTdy, 'k', 'LineWidth', 1.2);

% C) Dibujar el contorno de la placa
plot([0 4], [0 0.5], 'k-', 'LineWidth', 2); % Abajo
plot([0 4], [2 1.5], 'k-', 'LineWidth', 2); % Arriba
plot([0 0], [0 2], 'k-', 'LineWidth', 2); % Izquierda
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); % Derecha

% 7. Configuración Final
axis equal; % CRUCIAL para ver la ortogonalidad
axis([-0.5 4.5 -0.5 2.5]);
grid on;
xlabel('x'); ylabel('y');
title({'Isotermas y Gradiente \nablaT', ...
'T(x,y) = (1 + (y-1)^2)(4-x)'});
colorbar;

% --- CÁLCULO DEL MÁXIMO ---
% Busamos el valor máximo ignorando los NaNs
[max_temp, indice] = max(T(sonrisa));
[fil, col] = ind2sub(size(T), indice);
x_max = X(fil, col);
y_max = Y(fil, col);

fprintf('La temperatura máxima es %.2f en el punto (x=%.2f, y=%.2f)\n', ...
max_temp, x_max, y_max);
% Marcamos el punto máximo con una estrella roja
plot(x_max, y_max, 'p', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 15);
hold off;

4 Ley de Fourier

% 1. Configuración del Mallado
% Usamos h=0.2 para que las flechas no salgan amontonadas y se vea limpio
h = 0.2;
x = 0 : h : 4;
y = 0 : h : 2;
[X, Y] = meshgrid(x, y);

% 2. Definir la Geometría (La placa trapezoidal)
y_abajo = X ./ 8;
y_arriba = 2 - (X ./ 8);

% Máscara lógica: determina qué puntos están DENTRO de la placa
dentro = (Y >= y_abajo) & (Y <= y_arriba);

% 3. Definir la Temperatura T(x,y)
% Fórmula dada: T(x,y) = (1 + (y-1)^2)(4-x)
T_func = @(x,y) (1 + (y - 1).^2) .* (4 - x);
T = T_func(X, Y);

% 4. Cálculos Físicos (Ley de Fourier)
% A) Calculamos el Gradiente (nabla T) mediante derivadas numéricas
[dTdx, dTdy] = gradient(T, h, h);

% B) Calculamos el Flujo de Calor Q
% Q = -k * nabla T. (Suponiendo k=1) -> Q apunta al frío
Qx = -dTdx;
Qy = -dTdy;

% 5. Limpieza de datos (Recortar la figura)
% Asignamos NaN (Not a Number) a todo lo que esté fuera de la placa.
% Esto hace que MATLAB no dibuje nada en esas zonas.
Qx(~dentro) = NaN;
Qy(~dentro) = NaN;
T(~dentro) = NaN;

% 6. Visualización
figure('Color', 'w', 'Name', 'Flujo de Calor Fourier');
hold on;

% A) FONDO: Mapa de Calor (contourf)
% 'LineStyle','none' quita las rayitas negras entre colores para que sea suave.
contourf(X, Y, T, 20, 'LineStyle', 'none');
colormap('jet'); % Paleta de colores: Azul (frío) -> Rojo (caliente)

% Barra de color explicativa
c = colorbar;
c.Label.String = 'Temperatura T(x,y)';
c.Label.FontSize = 10;

% B) PROTAGONISTA: El Campo Vectorial Q (Flechas)
% Las pintamos de negro ('k') para que resalten sobre el fondo de color.
quiver(X, Y, Qx, Qy, 'k', 'LineWidth', 1.2, 'AutoScaleFactor', 1.5);

% C) DIBUJO DEL CONTORNO (Bordes Negros)
plot([0 4], [0 0.5], 'k-', 'LineWidth', 2); % Borde inferior
plot([0 4], [2 1.5], 'k-', 'LineWidth', 2); % Borde superior
plot([0 0], [0 2], 'k-', 'LineWidth', 2); % Borde izquierdo
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); % Borde derecho

% 7. Configuración Final
axis equal; % OBLIGATORIO: Para que los ángulos sean reales (ortogonalidad)
axis([-0.5 4.5 -0.5 2.5]); % Márgenes de visualización
grid on;
xlabel('Eje x'); ylabel('Eje y');
title({'Vector Flujo de Calor', ...
'El calor fluye desde las capas más calientes hacia las más frías'});

hold off;

5 Punto de máxima temperatura

% 1. Definición del Mallado y la Región
h = 0.05; % Paso del mallado (puedes ajustarlo a h=0.1 o h=0.01)
x = 0:h:4;
y = 0:h:2; % El rango de Y va de 0 a 2 según las funciones f(x) y g(x)
[X, Y] = meshgrid(x, y);

% Definición de los límites de la placa (Enunciado J)
% f(x) = x/8;
% g(x) = 2 - x/8;
Condition = (Y >= X./8) & (Y <= 2 - X./8);

% Aplicamos la máscara (Ponemos NaN fuera de la placa)
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo de las Derivadas Parciales (Gradiente)
% Usamos las fórmulas analíticas derivadas de T(x,y)
DerivadaX = -(1 + (Y - 1).^2);
DerivadaY = 2 .* (4 - X) .* (Y - 1);

% 3. Cálculo del Módulo del Gradiente
GradienteModulo = sqrt(DerivadaX.^2 + DerivadaY.^2);

% 4. Encontrar el Máximo Global
% Buscamos el valor máximo ignorando los NaNs
[MaxGradiente, IndiceLineal] = max(GradienteModulo);

% Convertimos el índice lineal a coordenadas de matriz (fila, columna)
[fila, col] = ind2sub(size(GradienteModulo), IndiceLineal);

% Obtenemos las coordenadas (x, y) y las componentes del vector en ese punto
PuntoMax = [X(fila, col), Y(fila, col)];
DirMax = [DerivadaX(fila, col), DerivadaY(fila, col)];

% Mostramos los resultados en la consola
fprintf('La variación máxima es: %.4f\n', MaxGradiente);
fprintf('Ocurre en el punto: (%.2f, %.2f)\n', PuntoMax(1), PuntoMax(2));
fprintf('Dirección del gradiente (vector): [%.2f, %.2f]\n', DirMax(1), DirMax(2));

% 5. Gráfica (Señalar con punto rojo)
figure;
hold on;
axis equal;
grid on;
title('Punto de Máxima Variación de Temperatura');
xlabel('Eje X'); ylabel('Eje Y');
axis([-0.5 4.5 -0.5 2.5]); % Ajuste de visualización

% Dibujar el contorno de la placa (opcional, para referencia visual)
x_borde = 0:0.1:4;
plot(x_borde, x_borde/8, 'k', 'LineWidth', 1.5); % f(x)
plot(x_borde, 2 - x_borde/8, 'k', 'LineWidth', 1.5); % g(x)
line([0 0], [0 2], 'Color', 'k', 'LineWidth', 1.5); % Pared izquierda
line([4 4], [0.5 1.5], 'Color', 'k', 'LineWidth', 1.5); % Borde derecho

% A) Dibujar el punto rojo
plot3(PuntoMax(1), PuntoMax(2), 10, 'ro', 'MarkerSize', 10, 'LineWidth', 3, 'MarkerFaceColor', 'r');

% B) Dibujar el vector de dirección en ese punto (Quiver)
% Normalizamos el vector para que se vea bien en la gráfica
escala = 0.5;
quiver(PuntoMax(1), PuntoMax(2), DirMax(1), DirMax(2), ...
escala, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5);

text(PuntoMax(1)+0.2, PuntoMax(2), sprintf('Max: %.2f', MaxGradiente));

hold off;

6 Campo de desplazamientos

%% Apartados 5 y 6: Campo de Desplazamientos y Deformación
clc; clear; close all;

% ---------------------------------------------------------
% 1. Definición del Mallado (Igual que antes)
% ---------------------------------------------------------
h = 0.1; % Paso un poco más fino para que la malla se vea bonita
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la región (Placa J)
Condition = (Y >= X./8) & (Y <= 2 - X./8);

% Aplicamos la máscara (NaN fuera de la placa)
X(~Condition) = NaN;
Y(~Condition) = NaN;

% ---------------------------------------------------------
% 2. Cálculo del Campo de Desplazamientos (u_x, u_y)
% ---------------------------------------------------------
% Primero necesitamos las coordenadas polares de cada punto (rho, theta)
Rho = sqrt(X.^2 + Y.^2);
Theta = atan2(Y, X);

% Definimos la componente escalar u_theta dada por el enunciado
% u(rho, theta) = -1/20 * rho^2 * cos(theta) * e_theta
U_theta_val = -1/20 .* Rho.^2 .* cos(Theta);

% Convertimos el vector de base polar a base cartesiana
% u_x = -u_theta * sin(theta)
% u_y = u_theta * cos(theta)
Ux = -U_theta_val .* sin(Theta);
Uy = U_theta_val .* cos(Theta);

% ---------------------------------------------------------
% 3. Calcular la posición deformada
% ---------------------------------------------------------
% Posición final = Posición inicial + Desplazamiento
X_final = X + Ux;
Y_final = Y + Uy;

% ---------------------------------------------------------
% 4. Gráfica Comparativa (Subplot como en tu ejemplo)
% ---------------------------------------------------------
figure('Name', 'Deformación de la Placa', 'Color', 'w');

% --- PLACA ORIGINAL ---
subplot(1, 2, 1);
% Usamos mesh con Z=0 para dibujar la 'rejilla' plana
mesh(X, Y, zeros(size(X)), 'EdgeColor', 'b');
view(2); % Vista superior (2D)
axis equal; grid on;
title('Placa antes del desplazamiento');
xlabel('Eje X'); ylabel('Eje Y');
axis([-0.5 4.5 -0.5 2.5]);

% --- PLACA DEFORMADA ---
subplot(1, 2, 2);
mesh(X_final, Y_final, zeros(size(X)), 'EdgeColor', 'b');
view(2);
axis equal; grid on;
title('Placa después del desplazamiento');
xlabel('Eje X'); ylabel('Eje Y');
axis([-0.5 4.5 -0.5 2.5]);

% ---------------------------------------------------------
% Extra: Apartado 5 (Vectores de desplazamiento - Quiver)
% ---------------------------------------------------------
% Si necesitas mostrar las flechitas de movimiento en la placa original:
figure('Name', 'Campo Vectorial de Desplazamientos');
hold on; axis equal; grid on;
title('Campo de Desplazamientos');

% Dibujamos contorno original de referencia
plot(X(:), Y(:), '.k', 'MarkerSize', 1);

% Dibujamos los vectores (quiver)
% Ajustamos 'scale' para que las flechas no sean gigantes visualmente
quiver(X, Y, Ux, Uy, 0, 'r', 'LineWidth', 1.5);

xlabel('Eje X'); ylabel('Eje Y');
xlim([-0.5 4.5]); ylim([-0.5 2.5]);

7 Desplazamiento del sólido

% 1. Definición de Geometría y Desplazamientos
h = 0.1; % Paso del mallado
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la Región (Placa)
% La placa existe si Y >= X/8 Y Y <= 2 - X/8
Condition = (Y >= X./8) & (Y <= 2 - X./8);

% Aplicamos máscara (NaN fuera de la placa para no dibujar ahí)
X(~Condition) = NaN;
Y(~Condition) = NaN;

% Cálculo de Desplazamientos
% Coordenadas polares necesarias para la fórmula
Rho = sqrt(X.^2 + Y.^2);
Theta = atan2(Y, X);

% Valor escalar de u_theta
U_theta_val = -1/20 .* Rho.^2 .* cos(Theta);

% Proyección a cartesianas (Crucial: u_x = -u_theta*sin, u_y = u_theta*cos)
Ux = -U_theta_val .* sin(Theta);
Uy = U_theta_val .* cos(Theta);

% --- Posición Final Deformada ---
X_final = X + Ux;
Y_final = Y + Uy;

% 2. Gráfica Comparativa de Sólidos Mallados (Como en tu ejemplo)
figure('Name', 'Comparación Antes/Después', 'Color', 'w', 'Position', [100, 100, 1000, 500]);

% SUBPLOT 1: Placa Original
subplot(1, 2, 1);
hold on;
% Usamos mesh. Z=zeros(size(X)) es para que sea un plano 2D.
% 'EdgeColor', 'c' pone las líneas en color cian
mesh(X, Y, zeros(size(X)), 'EdgeColor', [0 0.7 0.9], 'FaceColor', 'none');

view(2); % Fuerza la vista superior 2D
axis equal; % Mantiene la proporción real de los ejes
grid on;
box on; % Añade el marco completo

title('Placa antes del desplazamiento', 'FontWeight', 'bold');
xlabel('Eje X');
ylabel('Eje Y');
axis([-1.5 4.5 0 2.5]); % Ajuste de ejes similar a tu ejemplo

% SUBPLOT 2: Placa Deformada
subplot(1, 2, 2);
hold on;
% Dibujamos la malla usando las coordenadas finales (X_final, Y_final)
mesh(X_final, Y_final, zeros(size(X)), 'EdgeColor', [0 0.7 0.9], 'FaceColor', 'none');

view(2);
axis equal;
grid on;
box on;

title('Placa después del desplazamiento', 'FontWeight', 'bold');
xlabel('Eje X');
ylabel('Eje Y');
axis([-1.5 4.5 0 2.5]); % Mismos ejes para poder comparar visualmente

sgtitle('Visualización del Sólido Mallado bajo Deformación'); % Título

8 Divergencia

%% Apartado 7: Divergencia del Campo de Desplazamientos
clc; clear; close all;

% 1. Definición del Mallado y Región
h = 0.05;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa (Enunciado J)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo Analítico de la Divergencia
% Como demostramos analíticamente: Div(u) = y / 20
Divergencia = Y ./ 20;

% 3. Encontrar Máximos, Mínimos y Nulos analíticamente en el mallado
% Buscamos el valor máximo y mínimo en la matriz ignorando NaNs
max_val = max(Divergencia(:));
min_val = min(Divergencia(:));

% Encontramos los índices donde ocurren
[f_max, c_max] = find(Divergencia == max_val);
[f_min, c_min] = find(Divergencia == min_val);

% Puntos (x,y) correspondientes
PuntosMax = [X(f_max), Y(f_max)];
PuntosMin = [X(f_min), Y(f_min)];

% 4. Mostrar resultados en consola
fprintf('--- Resultados del Análisis de Divergencia ---\n');
fprintf('Divergencia Máxima: %.4f (Ocurre en y máximas)\n', max_val);
fprintf('Divergencia Mínima: %.4f (Ocurre en y = 0)\n', min_val);
fprintf('La divergencia es nula solo donde y = 0.\n');

% 5. Gráfica de la Divergencia
figure('Name', 'Divergencia del Desplazamiento', 'Color', 'w');
hold on;

% Dibujamos el mapa de colores (contourf rellena de color)
% Usamos 20 niveles para que se vea suave el gradiente
contourf(X, Y, Divergencia, 20, 'LineColor', 'none');
colormap(jet); % Mapa de colores clásico (Azul=Bajo, Rojo=Alto)
colorbar; % Barra lateral para saber qué valor es cada color

% Dibujar contorno de la placa para que quede bonito
x_borde = 0:0.1:4;
plot(x_borde, x_borde/8, 'k', 'LineWidth', 2);
plot(x_borde, 2 - x_borde/8, 'k', 'LineWidth', 2);
line([0 0], [0 2], 'Color', 'k', 'LineWidth', 2);

% Señalar Máximo y Mínimo
plot(PuntosMax(:,1), PuntosMax(:,2), 'rp', 'MarkerSize', 12, 'MarkerFaceColor', 'r'); % Estrella roja
plot(PuntosMin(:,1), PuntosMin(:,2), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b'); % Círculo azul

% Etiquetas y títulos
title('Cambio de Volumen Local');
xlabel('Eje X'); ylabel('Eje Y');
axis equal;
axis([-0.5 4.5 -0.5 2.5]);

% Leyenda manual para entender los puntos
legend('Divergencia (Color)', 'Bordes', '', '', 'Máxima Expansión', 'Mínima/Nula', ...
'Location', 'eastoutside');

hold off;

9 Rotacional

% 1. Definición del Mallado y Región
h = 0.05;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 1. Definición del Mallado y Región
h = 0.05;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa (Trapecio)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo del Rotacional (Magnitud)
%

10 Tensor deformaciones

%% Apartado 9 (Versión 2D): Tensiones Normales
clc; clear; close all;

% 1. Definición de Geometría y Material
h = 0.1;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Placa J (Trapecio)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN; Y(~Condition) = NaN;

lambda = 1; mu = 1;

% 2. Cálculo de Componentes
% Deformaciones y Divergencia
Div_u = Y ./ 20;
Eps_xx = Y ./ 20;
Eps_yy = zeros(size(X));
% Eps_zz es 0

% Tensiones Normales (Sigma = lambda*Div*I + 2*mu*Eps)
Sig_xx = lambda .* Div_u + 2 .* mu .* Eps_xx; % 3y/20
Sig_yy = lambda .* Div_u + 2 .* mu .* Eps_yy; % y/20
Sig_zz = lambda .* Div_u; % y/20 (Por efecto Poisson)
Sig_xy = 2 .* mu .* (-X./40); % -x/20 (Cizalladura)

% 3. Gráfica 2D (Subplot)
figure('Name', 'Apartado 9: Tensiones Normales (2D)', 'Color', 'w', 'Position', [100, 100, 1200, 400]);

% Sigma XX
subplot(1, 3, 1);
contourf(X, Y, Sig_xx, 20, 'LineColor', 'none');
axis equal; axis([-0.5 4.5 -0.5 2.5]); colorbar;
title('\sigma_{xx} (Eje i)'); xlabel('X'); ylabel('Y');

% Sigma YY
subplot(1, 3, 2);
contourf(X, Y, Sig_yy, 20, 'LineColor', 'none');
axis equal; axis([-0.5 4.5 -0.5 2.5]); colorbar;
title('\sigma_{yy} (Eje j)'); xlabel('X'); ylabel('Y');

% Sigma ZZ
subplot(1, 3, 3);
contourf(X, Y, Sig_zz, 20, 'LineColor', 'none');
axis equal; axis([-0.5 4.5 -0.5 2.5]); colorbar;
title('\sigma_{zz} (Eje k)'); xlabel('X'); ylabel('Y');

sgtitle('Componentes Normales del Tensor de Tensiones');

11 Tensiones tangenciales

% 1. Definición del Mallado
h = 0.15; % Un paso un poco más grande para que las flechas no se amontonen
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa 
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo de la Tensión Tangencial (Sigma_xy)
% Recordamos del apdo 9: Sig_xy = 2*mu*eps_xy = 2*1*(-x/40) = -x/20
Sig_xy = -X ./ 20;

% 3. Definición del Vector Tangencial
% En la cara vertical (normal i), la tangencial es vertical (j).
% El vector es: (0, Sigma_xy)
U_vec = zeros(size(X)); % Componente X es 0
V_vec = Sig_xy; % Componente Y es el valor de la tensión

% 4. Gráfica
figure('Name', 'Apartado 10: Vector Tensión Tangencial', 'Color', 'w');

% Usamos quiver para dibujar las flechitas azules
% quiver(x, y, u, v)
quiver(X, Y, U_vec, V_vec, 'b', 'LineWidth', 1);

% Configuración de la vista
view(2); % Vista 2D superior
axis equal; % Ejes proporcionales
grid on; % Cuadrícula de fondo

% Títulos y Etiquetas
title('Tensiones tangenciales (Plano ortogonal a i)');
xlabel('Eje X');
ylabel('Eje Y');

% Ajuste de límites para ver bien la placa
axis([-0.5 4.5 -0.5 2.5]);

12 Tensión de Von Mises

%% Apartado 11: Tensión de Von Mises (Usando autovalores 'eig')
% Inicializamos la matriz de resultados con NaNs
VonMises_Calc = NaN(size(X));
PuntoMax_VM = [-1, -1];
Max_VM_Val = -1;

% Bucle para recorrer cada punto del mallado
[filas, cols] = size(X);
for f = 1:filas
for c = 1:cols
% Si el punto está dentro de la placa (no es NaN)
if ~isnan(X(f,c))
% 1. Construir el tensor de tensiones en ese punto (3x3)
s_xx = Sig_xx(f,c);
s_yy = Sig_yy(f,c);
s_zz = Sig_zz(f,c);
s_xy = Sig_xy(f,c);

Tensor = [s_xx, s_xy, 0;
s_xy, s_yy, 0;
0, 0, s_zz];

% 2. Calcular autovalores (Tensiones Principales)
sig_principales = eig(Tensor);
s1 = sig_principales(1);
s2 = sig_principales(2);
s3 = sig_principales(3);

% 3. Fórmula de Von Mises dada en el enunciado
vm = sqrt( ((s1-s2)^2 + (s2-s3)^2 + (s3-s1)^2) / 2 );
VonMises_Calc(f,c) = vm;

% Guardar máximo para señalarlo después
if vm > Max_VM_Val
Max_VM_Val = vm;
PuntoMax_VM = [X(f,c), Y(f,c)];
end
end
end
end

% Gráfica
figure('Name', 'Apartado 11: Von Mises con Autovalores', 'Color', 'w');
contourf(X, Y, VonMises_Calc, 30, 'LineColor', 'none');
colormap(hot); % 'hot' suele usarse para mapas de "calor/tensión"
c = colorbar; c.Label.String = 'Tensión Von Mises (Pa)';

hold on;
axis equal; axis([-0.5 4.5 -0.5 2.5]);
title('Tensión de Von Mises (Cálculo mediante Autovalores)');
xlabel('Eje X'); ylabel('Eje Y');

% Señalar el punto máximo
plot(PuntoMax_VM(1), PuntoMax_VM(2), 'bp', 'MarkerSize', 14, 'MarkerFaceColor', 'b'); % Estrella azul
text(PuntoMax_VM(1)-1, PuntoMax_VM(2), sprintf('Máx: %.4f en (%.1f, %.1f)', Max_VM_Val, PuntoMax_VM), 'Color', 'b', 'FontWeight', 'bold');
hold off;

13 Campo de Fuerzas

%% Apartado 12: Campo de Fuerzas Volumétricas (F = -div Sigma)
clc; clear; close all;

% 1. Definición del Mallado
h = 0.1;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa (Trapecio)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

lambda = 1; mu = 1;

% 2. Recalculamos las Tensiones (Necesarias para derivar)
% Usamos las fórmulas analíticas que ya conocemos
Div_u = Y ./ 20;
Eps_xx = Y ./ 20;
Eps_yy = zeros(size(X));

Sig_xx = lambda .* Div_u + 2 .* mu .* Eps_xx; % = 3y/20
Sig_yy = lambda .* Div_u + 2 .* mu .* Eps_yy; % = y/20
Sig_xy = -X ./ 20; % = -x/20 (cizalladura)

% 3. Cálculo de la Divergencia del Tensor (Numéricamente)
% F = -div(sigma)
% Fx = -( d(Sxx)/dx + d(Sxy)/dy )
% Fy = -( d(Sxy)/dx + d(Syy)/dy )

% La función gradient calcula derivadas centrales.
% Sintaxis: gradient(Matriz, paso_x, paso_y)
[dSxx_dx, dSxx_dy] = gradient(Sig_xx, h, h);
[dSxy_dx, dSxy_dy] = gradient(Sig_xy, h, h);
[dSyy_dx, dSyy_dy] = gradient(Sig_yy, h, h);

Fx = - (dSxx_dx + dSxy_dy);
Fy = - (dSxy_dx + dSyy_dy);

% 4. Comprobación de Magnitud (Para ver si es cero)
MagnitudF = sqrt(Fx.^2 + Fy.^2);
max_F = max(MagnitudF(:));

fprintf('--- Resultado Apartado 12 ---\n');
fprintf('La fuerza máxima calculada es: %e\n', max_F);

if max_F < 1e-10
fprintf('CONCLUSIÓN: La fuerza es prácticamente NULA (Equilibrio).\n');
subtitle_text = 'Fuerza \approx 0 (Equilibrio Estático)';
else
subtitle_text = 'Existen fuerzas volumétricas';
end

% 5. Gráfica
figure('Name', 'Apartado 12: Campo de Fuerzas', 'Color', 'w');
hold on;

% Dibujamos los vectores.
% Como son casi cero, MATLAB pintará puntos. Para ver "algo",
% a veces se normalizan, pero lo correcto es mostrar que no hay fuerza.
quiver(X, Y, Fx, Fy, 'r');

% Dibujar bordes
plot([0 4 4 0 0], [0 0.5 1.5 2 0], 'k--', 'LineWidth', 0.5); % Caja aproximada

title('Campo de Fuerzas Volumétricas');
subtitle(subtitle_text);
xlabel('Eje X'); ylabel('Eje Y');
axis equal; axis([-0.5 4.5 -0.5 2.5]);
grid on;

% Añadimos una nota en la gráfica explicándolo
text(1, 1, {'Fuerza Nula', '(puntos indican magnitud ~0)'}, ...
'Color', 'r', 'FontWeight', 'bold', 'HorizontalAlignment', 'center');

hold off;

14 Densidad

%% Apartado 13: Cálculo de la Masa Total
clc; clear; close all;

% 1. Definición del Mallado
% Usamos un paso fino (h=0.01) para que la integral sea precisa
h = 0.01;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% 2. Definición de la Placa (Región)
Condition = (Y >= X./8) & (Y <= 2 - X./8);

% 3. Definición de la Densidad
% Función dada: d(x,y) = (4 - x) *