Estudio de la temperatura y deformación sobre una presa triangular

Trabajo realizado por estudiantes
Título	Estudio de la temperatura y deformación sobre una presa triangular
Asignatura	Teoría de Campos
Curso	2024-25
Autores	Jose Lema Mirella Espinal Arias Jorge García Sara López Caro Iria Cobos
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Se considera la sección transversal de una presa triangular, estando la superficie dada en coordenadas cartesianas, donde (x,y) pertenecen al siguiente recinto :

[math] [0,2]X[0,f(x)] [/math].

Siendo f(x) igual a: [math] f(x)=min(3,\frac{3}{2}(2-x)) [/math]

El objetivo de este estudio es analizar las diferentes transformaciones y aplicaciones que se le pueden ejercer a dicha presa triangular. Para ello se tienen que tener en cuenta tres funciones, la de temperatura, el campo de los desplazamientos y la función densidad de la presa.

• La temperatura viene dada por la función:

[math]T(x,y)=\frac{yx^2}{2}[/math]

• Los desplazamientos se corresponden con el campo:

[math] \vec u(x,y)=\frac{2(2-x)y\vec i-y\vec j}{50}[/math]

• Tomar como densidad:

[math]d(x, y) = (2 − |x − \frac{1}{2}|)(4 − y)[/math]

1 Mallado del solido

En primer lugar se tendrá que visualizar el mallado del solido. Para ello deberemos tener en cuenta el reciento planteado en el enunciado:

• Dibujar la presa triangular sabiendo que la superficie es [0,2] x [0,f(x)]
• Tener en cuenta los ejes de comando pedidos que en este caso son [-2,2] x [0,3]

Mallado de la Presa

h = 1/10; % Paso de Muestreo
x1 = 0:h:2;
y1 = 0:h:3;
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
Region = (y2 <= f(x2));
z2 = zeros(size(x2));
x2(~Region) = NaN; % Filtración de puntos fuera de la región(NaN)
y2(~Region) = NaN;
mesh(x2, y2, z2);
hold on;
x_borde = x1;
y_borde = arrayfun(f, x_borde);
plot3(x_borde, y_borde, zeros(size(x_borde)), 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
plot3([0 0], [0 f(0)], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2); %Delimitado de la superficie
plot3([2 2], [0 f(2)], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2); 
plot3([0 2], [0 0], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);  
axis equal;
axis([-1, 3, -1, 4]); %Asignacion del eje
xlabel('Eje X');
ylabel('Eje Y');
title('Mallado de la Presa');
view(2); 
grid on;
hold off;

2 Curvas de nivel y gradiente de temperatura

Las curvas de nivel son la representación de la temperatura sobre la presa, conociendo la función temperatura podemos calcular su dirección de crecimiento y punto máximo:

[math]T(x,y)=\frac{yx^2}{2}[/math].

2.1 Representación de temperatura y Punto Máximo

Temperatura Representada en Superficie

h = 1/10;%Paso muestreo
x1 = 0:h:2;        
y1 = 0:h:3;       
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2)*(2 - x));
region = (y2 <= f(x2)); 
Temp = (y2 .* x2.^2) / 2; %Funcion Temperatura
Temp(~region) = NaN; %Filtrar Sperficie Valida
surf(x2, y2, Temp, 'EdgeColor', 'none');   
hold on;
[maxTemp, idx] = max(Temp(:)); % Máximo valor 
[x_max, y_max] = ind2sub(size(Temp), idx); 
x_coord = x2(x_max, y_max); % Coordenada X del máximo
y_coord = y2(x_max, y_max); % Coordenada Y del máximo
plot3(x_coord, y_coord, maxTemp, 'ro', 'MarkerSize', 8, 'LineWidth', 2); % Punto rojo en 3D de coordenada maxima
text(x_coord, y_coord, maxTemp, sprintf('  Maximo: %.2f', maxTemp), 'Color', 'r', 'FontSize', 10);
axis equal;
axis([0, 2, 0, 3]);  
xlabel('Eje X');
ylabel('Eje Y');
title('Temperatura representada sobre la Superficie');
colorbar;  
grid on;  
hold off;

2.2 Representacion Curvas de nivel de T

Curvas de Nivel de la Presa

h = 1/10;
x1 = 0:h:2;       
y1 = 0:h:3;        
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2)*(2 - x));
region = (y2 <= f(x2)); 
T = (y2 .* x2.^2) / 2;
T(~region) = NaN;
contour3(x2, y2, T, 20, 'LineWidth', 2); %Mostrar Curvas de Nivel
colorbar;                                 
axis equal;
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Temperatura T(x, y)');
title('Curvas de Nivel de la Temperatura representadas');
grid on;
hold off;

2.3 Gradiente de T

El gradiente de la temperatura \(T(x,y)\) define en que dirección el campo de temperaturas varia mas rápido, este se define de la siguiente manera:

[math]\nabla T(x,y) =\frac{∂T}{∂x}\vec i + \frac{∂T}{∂y}\vec j + \frac{∂T}{∂z}\vec ez [/math]

Conociendo dicha definición en este caso el gradiente será:

[math]\nabla T(x,y) = yx\vec i + \frac{x^2}{2}\vec j [/math]

h = 1/10;
x1 = 0:h:2;        
y1 = 0:h:3;       
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
region = (y2 <= f(x2)); 
T = (y2 .* x2.^2) / 2;
T(~region) = NaN;
[Tx, Ty] = gradient(T, h);%Gradiente de la Temperatura
Tx(~region) = NaN;%Filtrar Gradiente fuera de la Región
Ty(~region) = NaN;
z_surface = T;
contour3(x2, y2, z_surface, 20, 'LineWidth', 2);
hold on;
scale = 3; 
quiver3(x2, y2, z_surface, scale*Tx, scale*Ty, zeros(size(Tx)), 'r', 'LineWidth', 1.5);
axis equal;
axis([0, 2, 0, 3, 0, 1]); 
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Temperatura T(x, y)');
title('Curvas de Nivel y Gradiente de la Temperatura');
colorbar;
grid on;
hold off;

3 Ley de Fourier

La ley de Fourier establece que la energía calorífica [math]\vec Q~[/math] viaja a través de la superficie con la siguiente formula:

[math]\vec Q=-k∇T[/math]

Sabiendo que K es la constante de conductividad térmica de la presa que supondremos k=1:

[math]\vec Q=-k∇T=-yx\vec i -\frac{x^2}{2}\vec j[/math]

h = 1/10;
x1 = 0:h:2;        
y1 = 0:h:3;        
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
region = (y2 <= f(x2)); 
T = (y2 .* x2.^2) / 2;
T(~region) = NaN;
[Tx, Ty] = gradient(T, h);
Qx = -Tx;%Flujo de la energía calorifica
Qy = -Ty;
Qx(~region) = NaN;%Filtrar flujo en superficie
Qy(~region) = NaN;
z_surface = T;
contour3(x2, y2, z_surface, 20, 'LineWidth', 2);
hold on;
scale = 1.5; 
quiver3(x2, y2, z_surface, scale*Qx, scale*Qy, zeros(size(Qx)), 'r', 'LineWidth', 1.5);
axis equal;
axis([0, 2, 0, 3, 0, 1]); 
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Temperatura T(x, y)');
title('Curvas de Nivel y Flujo de Energía Calorífica en 3D');
colorbar; 
grid on;
hold off;

4 Máxima temperatura

Es el punto en el dominio donde el valor de la función [math]T(x,y)=\frac{yx^2}{2}[/math] es máximo. La temperatura máxima es 0.8820ºC.

Máxima temperatura

h = 1/10;  
x1 = 0:h:2;
y1 = 0:h:3;
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
Region = (y2 <= f(x2));  
z2 = zeros(size(x2));
x2(~Region) = NaN;
y2(~Region) = NaN;
T = (y2 .* x2.^2) / 2; % Función de temperatura T(x, y)
% Gradiente de T(x, y)
dTdx = y2 .* x2;       % Derivada parcial respecto a x
dTdy = x2.^2 / 2;      % Derivada parcial respecto a y
grad = sqrt(dTdx.^2 + dTdy.^2);  
[max_grad, idx] = max(grad_magnitude(:)); % Punto de mayor gradiente
[max_x, max_y] = ind2sub(size(grad_magnitude), idx);
max_point = [x2(max_x, max_y), y2(max_x, max_y)];
figure;
hold on;
mesh(x2, y2, z2);
x_borde = x1;
y_borde = arrayfun(f, x_borde);
plot3(x_borde, y_borde, zeros(size(x_borde)), 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2); 
plot3([0 0], [0 f(0)], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
plot3([2 2], [0 f(2)], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
plot3([0 2], [0 0], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
% Representar direcciones de gradiente
quiver3(x2, y2, z2, dTdx, dTdy, zeros(size(dTdx)), 'k', 'LineWidth', 1.5);
% Resaltar el punto de mayor gradiente
plot3(max_point(1), max_point(2), 0, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
axis equal;
axis([-1, 3, -1, 4, -1, 1]);
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
title('Mallado de la Superficie y Dirección de Máxima Variación de Temperatura');
view(3); % Vista 3D
grid on;
hold off;
Tempmax=max(max(T))

5 Campo de vectores

En este apartado mostraremos el campo de vectores en los puntos del mallado del solido. Nuestro campo es: [math] \vec u(x,y)=\frac{2(2-x)y\vec i-y\vec j}{50}[/math].

En el campo de desplazamientos que se observa en la imagen, los puntos que están fijos son aquellos donde el vector de desplazamiento es cero.

Campo de vectores

h = 1/10;  
x1 = 0:h:2;
y1 = 0:h:3;
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
Region = (y2 <= f(x2));  
z2 = zeros(size(x2));
x2(~Region) = NaN; 
y2(~Region) = NaN;
% Campo de desplazamientos u(x, y)
ux = (2 * (2 - x2) .* y2) / 50;  % Componente x
uy = -y2 / 50;                  % Componente y
scale = 1.5;
figure;
hold on;
mesh(x2, y2, z2, 'EdgeColor', [0.8 0.8 0.8], 'FaceAlpha', 0.2); 
quiver3(x2, y2, z2, ux, uy, zeros(size(ux)), scale, 'k'); % Flechas de desplazamiento
axis equal;
axis([-0.5, 2.5, -0.5, 3.5, -0.5, 0.5]);
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
title('Mallado de la Superficie y Campo de Desplazamientos');
view(3); % Vista 3D
grid on;
hold off;

6 Sólido antes y después del desplazamiento

Estudiamos el desplazamiento de una presa triangular al campo [math] \vec u(x,y)=\frac{2(2-x)y\vec i-y\vec j}{50}[/math]

6.1 Antes del desplazamiento

Campo de vectores

h = 1/10; 
x1 = 0:h:2; 
y1 = 0:h:3; 
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
Region = (y2 <= f(x2)); 
z2 = zeros(size(x2));
x2(~Region) = NaN;
y2(~Region) = NaN;
figure;
surf(x2, y2, z2, 'FaceColor', [0.2, 0.6, 0.8], 'EdgeColor', [0.2, 0.2, 0.2]); % Superficie con líneas internas visibles
hold on;
x_contorno = [0, 2, 0, 0];
y_contorno = [0, 0, 3, 0];
plot(x_contorno, y_contorno, 'b-', 'LineWidth', 2); % Línea del contorno en azul
axis equal;
axis([-1, 3, -1, 4]);
xlabel('Eje X');
ylabel('Eje Y');
title('Presa triangular');
view(2); 
grid on;
hold off;

6.2 Después del desplazamiento

Campo de vectores

h = 1/10;
x1 = 0:h:2; 
y1 = 0:h:3;
[x2, y2] = meshgrid(x1, y1);
f = @(x) min(3, (3/2) * (2 - x));
Region = (y2 <= f(x2)); 
z2 = zeros(size(x2));
x2(~Region) = NaN;
y2(~Region) = NaN;
% Campo de desplazamientos u(x, y)
ux = (2 * (2 - x2) .* y2) / 50; % Componente x del desplazamiento
uy = -y2 / 50;                  % Componente y del desplazamiento
% Aplicar el desplazamiento
x2_desplazado = x2 + ux;
y2_desplazado = y2 + uy;
figure;
surf(x2_desplazado, y2_desplazado, z2, ...
    'FaceColor', 'magenta', 'EdgeColor', [0.2, 0.2, 0.2], 'FaceAlpha', 0.7);
hold on;
x_contorno = [0, 2, 0, 0];
y_contorno = [0, 0, 3, 0];
% Aplicar el desplazamiento al contorno
ux_contorno = (2 * (2 - x_contorno) .* y_contorno) / 50;
uy_contorno = -y_contorno / 50;
x_contorno_desplazado = x_contorno + ux_contorno;
y_contorno_desplazado = y_contorno + uy_contorno;
plot(x_contorno_desplazado, y_contorno_desplazado, 'b-', 'LineWidth', 2);
axis equal;
axis([-1, 3, -1, 4]);
xlabel('Eje X');
ylabel('Eje Y');

7 Divergencia del Vector U

El campo vectorial desplazamiento [math] \vec u(x,y)=\frac{2(2-x)y\vec i-y\vec j}{50}[/math] es igual a [math] \vec u(x,y)=(\frac{4 y}{50} - \frac{2 x y}{50})\vec{i} - (\frac{y}{50})\vec{j} [/math] A partir de esa expresión se calcula la divergencia: [math]\ \nabla \cdot \vec u =\frac{\partial \vec u_1}{\partial x} + \frac{\partial \vec u_2}{\partial y} = \frac{\partial (\frac{4y}{50})}{\partial x} - \frac{\partial (\frac{2xy}{50})}{\partial x} + \frac{\partial(\frac{-y}{50})}{\partial y} = \frac{-2y -1}{50}[/math]

El máximo de la divergencia se encuentra en -0.02, donde la divergencia es menos negativa y el mínimo de la divergencia en -0.14 donde el campo vectorial tiene la máxima contracción en ese sistema.

La divergencia es nula en todos los puntos donde [math]\ \nabla \cdot \vec u=-\frac{1}{2}[/math] para cualquier valor de [math]\ x [/math].

Divergencia

h = 1/10; 
x1 = 0:h:2; 
y1 = 0:h:3; 
f = @(x) min(3, (3/2) * (2 - x)); 
[x2, y2] = meshgrid(x1, y1); 
Region = (y2 <= f(x2)); 
z2 = zeros(size(x2)); 
x2(~Region) = NaN;
y2(~Region) = NaN;

% Cálculo de la divergencia
D = (-2 .* y2 / 50) - (1 / 50);
% Graficar el mallado interno y divergencia
figure;
% Superficie de divergencia con mallado interno
mesh(x2, y2, D, 'FaceColor', 'interp', 'EdgeColor', [0.5, 0.5, 0.5]); 
hold on;
x_borde = x1;
y_borde = arrayfun(f, x_borde);
plot3(x_borde, y_borde, zeros(size(x_borde)), 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2); 
plot3([0 0], [0 f(0)], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
plot3([2 2], [0 f(2)], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
plot3([0 2], [0 0], [0 0], 'Color', [0.2, 0.4, 0.6], 'LineWidth', 2);
mesh(x2, y2, z2, 'EdgeColor', 'black', 'FaceColor', 'none'); 
axis equal;
axis([-1, 3, -1, 4, min(D(:)), max(D(:))]);
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Divergencia');
title('Divergencia');
colorbar; % Barra de colores para la divergencia
view(3); 
grid on;
hold off;
% Cálculo de máximos y mínimos
Dmax = max(D(:), [], 'omitnan'); % Máximo de la divergencia
Dmin = min(D(:), [], 'omitnan'); % Mínimo de la divergencia
disp(['Máximo de la divergencia: ', num2str(Dmax)]);
disp(['Mínimo de la divergencia: ', num2str(Dmin)]);

8 Rotacional

Para empezar se calcula el rotacional de [math]\vec{u}[/math]:

[math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{i} & \vec{j} &\vec{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y}\ & \frac{\partial }{\partial z}\\ u_{1} & u_{2} & u_{3} \end{vmatrix}=\begin{vmatrix}\vec{i} & \vec{j} &\vec{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y}\ & \frac{\partial }{\partial z}\\ (\frac{4 y}{50} - \frac{2 x y}{50}) & - (\frac{y}{50}) & 0 \end{vmatrix}=(\frac{2 x }{50}-\frac{4}{50})\vec{k}[/math]

Calculando el módulo:

[math]|∇ × \vec{u}|= \frac{-2(2-x)}{50} [/math]

Rotacional

h = 1/10; 
x = 0:h:2; 
f = @(x) min(3, 3/2 * (2 - x)); 
y = 0:h:max(f(x)); 

% Crear el mallado
[MX, MY] = meshgrid(x, y);
MY(MY > f(MX)) = NaN; 

% Calcular componentes del rotacional
Ri = 0 * MY; 
Rj = 0 * MY; 
Rk = (-(2 - MX)) ./ 25; 

% Graficar el rotacional
figure;
quiver3(MX, MY, zeros(size(MX)), Ri, Rj, Rk, 'k', 'LineWidth', 1.5);
hold on;
mesh(MX, MY, zeros(size(MX)), 'EdgeColor', [0.2, 0.4, 0.6], 'FaceColor', 'none');
hold off;
axis equal;
axis([-0.5, 2.5, -0.5, 3.5, -0.5, 1]);
title('Rotacional');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
view(3);
grid on;

9 Tensor de deformaciones

Para este apartado hay que definir dos nuevos conceptos. El primero es el tensor de deformaciones: [math]Ԑ(\vec{u}) = \frac{1}{2}(∇\vec{u} + ∇\vec{u}^t)[/math] que es la parte simétrica de [math]\vec{u}[/math].

El otro es el tensor de tensiones: [math]σ=λ∇·\vec{u}1 + 2μԐ[/math], donde [math]1[/math] es el tensor identidad en el conjunto de vectores libres del espacio [math]R^3[/math] y [math]λ[/math] y [math]µ[/math] son coeficientes de Lamé que dependen de las propiedades elásticas de cada material.

En este apartado se desea comprobar si existen tensiones ortogonales a la dirección de la placa, y si existen, representarlas. Para ello, se toman los valores de [math]λ=µ=1[/math].

Se empieza calculando el gradiente de [math]\vec{u}[/math], para ello se calcula el gradiente del campo vectorial: [math] ∇\vec{u}=\begin{pmatrix}\frac{\partial u_{1}}{\partial x} & \frac{\partial u_{1}}{\partial y} & \frac{\partial u_{1}}{\partial t} \\ \frac{\partial u_{2}}{\partial x} & \frac{\partial u_{2}}{\partial y} & \frac{\partial u_{2}}{\partial t}\\ \frac{\partial u_{3}}{\partial x} & \frac{\partial u_{3}}{\partial y} & \frac{\partial u_{3}}{\partial t}\end{pmatrix} = \begin{pmatrix} - \frac{2y}{50}\ & \frac{4}{50}-\frac{2x}{50}\ & 0\\0 & - \frac{1}{50}\ & 0\\ 0 & 0 & 0\end{pmatrix}[/math];

También se debe calcular la traspuesta del gradiente de [math]\vec{u} [/math] : [math] ∇\vec{u}^t= \begin{pmatrix} -\frac{2y}{50} & 0 & 0 \\ \frac{4}{50}-\frac{2x}{50} & -\frac{1}{50} & 0\\ 0 & 0 & 0\end{pmatrix} [/math].

Con todo ello se puede calcular el tensor de deformaciones [math]Ԑ(\vec{u}) = \frac{∇\vec{u} + ∇\vec{u}t}{2}=\begin{pmatrix}- \frac{y}{25}\ & \frac{2(2-x)}{100}\ & 0\\\frac{2(2-x)}{100} & - \frac{1}{50}\ & 0\\ 0 & 0 & 0\end{pmatrix}[/math]

Ya calculado la divergencia anteriormente [math]\ \nabla \cdot \vec u = \frac{-2y -1}{50}[/math] se podrá sustituir en la fórmula del tensor de tensiones:

[math]σ=λ∇·\vec{u}1 + 2μԐ[/math]= [math] 1 · (-\frac{2y}{50} -\frac{1}{50})[/math] \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} + [math] 2 · 1 [/math] \begin{pmatrix}- \frac{2y}{50}\ & \frac{2-x}{50}\ & 0\\\frac{2-x}{50} & - \frac{1}{50}\ & 0\\ 0 & 0 & 0\end{pmatrix} = \begin{pmatrix} -\frac{6y+1}{50} & \frac{4-2x}{50} & 0\\ \frac{2-x}{25} & - \frac{y}{25}-\frac{3}{50} & 0\\ 0 & 0 & - \frac{2y}{50}-\frac{1}{50}\end{pmatrix}

Una vez se ha obtenido el tensor de tensiones, se calculan las tensiones normales respecto a cada uno de los ejes.

[math]\vec{i}·σ·\vec{i}=\begin{pmatrix} 1 & 0 & 0 \end{pmatrix}\begin{pmatrix} -\frac{6y+1}{50} & \frac{4-2x}{50} & 0\\ \frac{4-2x}{50} & - \frac{3+2y}{50} & 0\\ 0 & 0 & - \frac{2y}{50}-\frac{1}{50}\end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}=\begin{pmatrix} 1 & 0 & 0 \end{pmatrix}\begin{pmatrix} -\frac{3y}{25}-\frac{1}{50} \\ \frac{2-x}{25} \\ 0 \end{pmatrix}= -\frac{6y+1}{50}[/math]

[math]\vec{j}·σ·\vec{j}=\begin{pmatrix} 0 & 1 & 0 \end{pmatrix}\begin{pmatrix} -\frac{6y+1}{50} & \frac{4-2x}{50} & 0\\ \frac{4-2x}{50} & - \frac{3+2y}{50} & 0\\ 0 & 0 & - \frac{2y}{50}-\frac{1}{50}\end{pmatrix} \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}=\begin{pmatrix} 0 & 1 & 0 \end{pmatrix}\begin{pmatrix} \frac{2-x}{25} \\ -\frac{y}{25}-\frac{3}{50} \\ 0 \end{pmatrix}=-\frac{3+2y}{50}[/math]

[math]\vec{k}·σ·\vec{k}=\begin{pmatrix} 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} -\frac{6y+1}{50} & \frac{4-2x}{50} & 0\\ \frac{4-2x}{50} & - \frac{3+2y}{50} & 0\\ 0 & 0 & - \frac{2y}{50}-\frac{1}{50}\end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}=\begin{pmatrix} 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} 0 \\ 0 \\ -\frac{y}{25}-\frac{1}{50}\end{pmatrix}= -\frac{y}{25}-\frac{1}{50}[/math]

h = 1/10; 
x = 0:h:2;
f = @(x) min(3, 3/2 * (2 - x)); 
y = 0:h:max(f(x)); 
[MX, MY] = meshgrid(x, y);
MY(MY > f(MX)) = NaN;

Tni = (-3 .* MY) ./ 25 - 1 / 50; 
Tnj = -MY ./ 25 - 3 / 50;      
Tnk = -MY ./ 25 - 1 / 50;       

figure;
set(gcf, 'Position', [100, 100, 1200, 600]); % Configurar el tamaño de la ventana de la figura

% Tensión normal en dirección i
subplot(1, 3, 1);
surf(MX, MY, Tni);
shading interp;
axis equal;
xlabel('X', 'FontSize', 14);
ylabel('Y', 'FontSize', 14);
zlabel('Z', 'FontSize', 14);
view(2);
colorbar;
title('Tensión normal en i', 'FontSize', 14);

% Tensión normal en dirección j
subplot(1, 3, 2);
surf(MX, MY, Tnj);
shading interp;
axis equal;
xlabel('X', 'FontSize', 14);
ylabel('Y', 'FontSize', 14);
zlabel('Z', 'FontSize', 14);
view(2);
colorbar;
title('Tensión normal en j', 'FontSize', 14);

% Tensión normal en dirección k
subplot(1, 3, 3);
surf(MX, MY, Tnk);
shading interp;
axis equal;
xlabel('X', 'FontSize', 14);
ylabel('Y', 'FontSize', 14);
zlabel('Z', 'FontSize', 14);
view(2);
colorbar;
title('Tensión normal en k', 'FontSize', 14);

10 Tensiones Tangenciales

En este apartado, se calcularán las tensiones tangenciales respecto al plano ortogonal a u, en este caso el plano ortogonal a i. En la Figura se dibujaran solo las que no sean nulas. Seguiremos la fórmula: [math] |σ·\vec{i}−(\vec{i}·σ·\vec{i})\vec{i}|[/math].

[math]σ·\vec{i}=\begin{pmatrix} - \frac{3y}{25}-\frac{1}{50} & \frac{2-x}{25} & 0\\ \frac{2-x}{25} & - \frac{y}{25}-\frac{3}{50} & 0\\ 0 & 0 & - \frac{y}{25}-\frac{1}{50}\end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}=\begin{pmatrix} - \frac{3y}{25}-\frac{1}{50} \\ \frac{2-x}{25} \\ 0 \end{pmatrix} [/math]

[math] (\vec{i}·σ·\vec{i})\vec{i}= -\frac{3y}{25}-\frac{1}{50} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}=\begin{pmatrix} -\frac{3y}{25}-\frac{1}{50} \\ 0 \\ 0 \end{pmatrix} [/math]

[math] |σ·\vec{i}−(\vec{i}·σ·\vec{i})\vec{i}| = \begin{pmatrix} - \frac{3y}{25}-\frac{1}{50} \\ \frac{2-x}{25} \\ 0 \end{pmatrix} - \begin{pmatrix} -\frac{3y}{25}-\frac{1}{50} \\ 0 \\ 0 \end{pmatrix} = |\frac{2-x}{25} \vec{j}| = \frac{2-x}{25} [/math]

figura

h = 1/10;                 
x = 0:h:2;                 
f = @(x1) min(3, 3/2 * (2 - x1)); 
y = 0:h:max(f(x1));

[X, Y] = meshgrid(x, y); 
Y(Y > f(X)) = NaN;
tang=(2-X)./25;%Función de la tensión tangencial
subplot(2,1,1)  %Representación en 2D
surf(X,Y,tang)
shading interp
axis equal
xlabel('X')
ylabel('Y')
zlabel('Z')
view(2)
colorbar
title('Tensión tangencial a la dirección del eje e i')
subplot(2,1,2)  %Representación en 3D
surf(X,Y,tang)
shading interp
axis equal
xlabel('X')
ylabel('Y')
zlabel('Z')
view(3)
colorbar
title('Tensión tangencial a la dirección del eje e i')

11 Tensiones de Von Mises

Para la tensión de Von Mises seguiremos la siguiente fórmula:

• [math] \sigma_V= \sqrt \frac{(\sigma _1-\sigma _2)^2+(\sigma _2-\sigma _3)^2+(\sigma _3-\sigma _1)^2}{2} [/math]

Donde σ1, σ2 y σ3 (también conocidos como tensiones principales) son los autovalores de σ (tensor de tensiones). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro).

Apartado 11

h = 1/10;
x = 0:h:2;
f = @(x) min(3, 3/2 * (2 - x)); 
y = 0:h:max(f(x));
[X, Y] = meshgrid(x, y);

% Inicializar la matriz de Von Mises
VonMises = zeros(size(X));

% Definición de las funciones para las tensiones
M11 = @(x, y) (-6 .* y - 1) ./ 50;
M12 = @(x, y) (2 - x) / 25;
M22 = @(x, y) (-2 .* y - 3) ./ 50;

% Iterar sobre el mallado para calcular Von Mises en cada punto
for i = 1:size(Y, 1)
    for j = 1:size(X, 2)
        % Construcción de la matriz de tensiones
        sigma = [M11(X(i, j), Y(i, j)), M12(X(i, j), Y(i, j)), 0;
                 M12(X(i, j), Y(i, j)), M22(X(i, j), Y(i, j)), 0;
                 0, 0, (-2 .* Y(i, j) - 1) ./ 50];

        % Cálculo de los autovalores de la matriz de tensiones
        eigvals = eig(sigma);

        % Cálculo de la tensión de Von Mises
        VonMises(i, j) = sqrt(((eigvals(1) - eigvals(2))^2 + ...
                               (eigvals(2) - eigvals(3))^2 + ...
                               (eigvals(3) - eigvals(1))^2) / 2);
    end
end

% Limitar la gráfica a los puntos dentro del triángulo definido por f(x)
mask = Y > f(X);
VonMises(mask) = NaN;

% Graficar los resultados en tres vistas diferentes
subplot(1, 3, 1);
surf(X, Y, VonMises);
view(2);
xlabel('Eje x');
ylabel('Eje y');
axis([-0.5, 2.5, -0.5, 3.5]);
colorbar;
title('Plano XOY');

subplot(1, 3, 2);
surf(X, Y, VonMises);
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
axis([-0.5, 2.5, -0.5, 3.5, -0.5, 0.5]);
colorbar;
title('Tensión de Von Mises en 3D');
axis vis3d;

subplot(1, 3, 3);
hold on;
surf(X, Y, VonMises);
view(0, 0);
xlabel('Eje x');
zlabel('Eje z');
axis([-0.5, 2.5, -0.5, 3.5, -0.5, 0.5]);
colorbar;
title('Plano XOZ');
MAX = max(VonMises, [], 'all');
txt = ['Máxima tensión: ' num2str(MAX)];
text(-0.5, -0.5, MAX * 0.8, txt);
hold off;

12 Campo de Fuerzas Causantes del Desplazamiento

El Campo de Fuerzas que actúan en la placa plana triangular, que son causantes de los desplazamiento, se aproximan usando la ecuación de la elasticidad lineal, modelo matemático que define como los solidos deformables responden a las fuerzas externas a las que son sometidas, manteniendo una relación lineal entre tensión y deformación. Con esta teoría se asume que el solido es isótropo y homogéneo, es decir que tiene las mismas características mecánicas en todos los puntos de su estructura, en este caso en todos los puntos de la estructura de la Presa Triangular.

Con esta formula definimos el Campo de Fuerzas que causan el desplazamiento: [math]\vec{F}=-∇·σ[/math]

Recordando σ del apartado 9, se calcula el campo de fuerzas pedido con el siguiente comando:

campodefuerzas

h = 1/10;  
x = 0:h:2; 
y = 0:h:3; 
[X, Y] = meshgrid(x, y);
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;  

F1 = (-3*Y)/25 - 1/25 + (4*(2 - X))/25;
F2 = (4*(2 - X))/25 + (-Y - 2)/25;
[div_x1, div_y_1] = gradient(F1, x, y);
[div_2x, div_2y] = gradient(F2, x, y);
divergencia = div_x1 + div_2y;
Fx = -divergencia;
Fy = -divergencia;
Fuerza = sqrt(Fx.^2 + Fy.^2);
max_F = max(Fuerza(:));
figure;
quiver(X, Y, Fx, Fy);
xlabel('x');
ylabel('y');
title('Fuerza en la presa triangular');

13 Calculo de la Masa

Definimos que la presa triangular tiene una densidad definida por: [math]d(x, y) = (2 − |x − \frac{1}{2}|)(4 − y)[/math]

Para calcular la masa cuando una densidad depende de 2 variables, como es el caso, integramos la densidad sobre la region correspondiente en el plano XY.

Definimos las regiones de la Presa Triangular:

Coordenadas X en [0;2]

Coordenadas Y en [0;f(x)]

Definimos la funcion F(x): [math] f(x)=min(3,\frac{3}{2}(2-x)) [/math]

Donde con ayuda de Matlab, se define que si integramos la densidad en las regiones por las que esta definida la presa Triangular.

[math] M = \int \int_{A} d(x, y) \, dx \, dy [/math]

Remplazamos en la ecuacion la Region a integrar y la funcion de Densidad Definida anteriormente.

[math] M = \int_{0}^{2} \int_{0}^{f(x)} \left( 2 - \left| x - \frac{1}{2} \right| \right)(4 - y) \, dy \, dx [/math]

Posteriormente se realiza la parametrizacion de la Region a integrar y hallamos el valor absoluto. [math]\vec{r'_u}=\vec{j}[/math]

[math]\vec{r'_v}=\vec{i}[/math]

[math]\vec{r'_u} × \vec{r'_v}[/math] = [math]\begin{vmatrix} \vec{i} & \vec{j} & \vec{k}\\ 1 & 0 & 0 \\ 0 & 1 & 0\end{vmatrix} = \vec{k}; [/math]

[math]|\vec{r'_u} × \vec{r'_v}|=1[/math]

Luego se hace reemplazan las variables y se aplica la definicion de la Integral para calcular la Masa.

[math]\int_{A} d( \vec{r} (u,v) |\vec{r'_u} × \vec{r'_v}| dudv=\int_{0}^{2}\int_{0}^{f(x)} d( \vec{r} (u,v) |\vec{r'_u} × \vec{r'_v}| dudv=\int_{0}^{2}\int_{0}^{f(x)} (2 − |u − \frac{1}{2}|)(4 − v) dudv= 14.19 [/math]

La masa Total es 14.19

% Definir la función límite superior f(x)
f = @(x) min(3, (3/2) * (2 - x)); % Función f(x)

% Definir la densidad d(x, y)
d = @(x, y) (2 - abs(x - 1/2)) .* (4 - y); % Densidad

% Calcular la masa usando integración doble
masa = integral(@(x) arrayfun(@(x) integral(@(y) d(x, y), 0, f(x)), x), 0, 2);

% Mostrar el resultado
disp(['La masa del objeto es: ', num2str(masa)]);