Campos de temperaturas y deformaciones en un cuarto de anillo (32)

Trabajo realizado por estudiantes
Título	Campos de temperaturas y deformaciones en una placa 2-d. Grupo 32
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Maleno Ayala, Rubén Sánchez Gil, Miguel González López, Mario Ruiz López, Álvaro
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Un Campo de Temperaturas es un Campo Escalar aplicado a un objeto físico con el objetivo de representar las temperaturas propias de los puntos del objeto.

La representación visual de un campo escalar puede hacerse mediante programas como Octave y MATLAB, los cuales utilizan operaciones escalares y elementos matriciales para realizar cálculos complejos.

1 Mallado de la placa

Para generar el mallado del cuarto de anillo circular usaremos el comando meshgrid en Matlab. El cuarto de anillo está centrado en \((0,0)\) y sus radios son \(R_{1} = 1\) y \(R_{2} = 2\). Además, se encuentra en el plano \(y ≥\frac{|x|}{2}\).

Su representación está encuadrado en la región \((x, y) ∈ [-3; 3] × [−1; 3]\). El paso de muestreo se refiere a las subdivisiones deseadas por unidad en ambos ejes, en este caso, \(h = 2/10\).

Ahora pasaremos estas características a coordenadas cilíndricas, donde \(ρ\) conservará el radio del anillo y \(θ\) dependerá del ángulo que forme las rectas \(y =\frac{x}{2}\) e \(y =\frac{-x}{2}\) con los ejes \(X\) e \(Y\), cuyo valor es \(α = tg(\frac{1}{2}) ≈ 26,6°\), por tanto:

[math](ρ,θ) ∈ [1; 2] × [α,π-α][/math]

Figura 1

clc
% Definición de variables
paso_muestreo = 2/10;
rho = 1:paso_muestreo:2;
alfa = atan(1/2);
theta = linspace(a,pi-alfa,pi/paso_muestreo);
[RHO,THETA] = meshgrid(rho,theta);
% Conversión a cilíndricas
x = RHO.*cos(THETA);
y = RHO.*sin(THETA);
mesh(x,y,0.*x)
view(2)
axis([-3,3,-1,3])
% Definición del título
title('CUARTO DE ANILLO CIRCULAR')

2 Curvas isotermas

Una vez hemos representado el sólido, podemos dibujar las curvas de nivel de la temperatura. La distribución de la temperatura en el sólido viene dada por el campo escalar: [math]T(x,y)=cos(x^2)+sin((y−1)^2)[/math]

A partir de este campo escalar, podemos definir el gradiente de T como [math]\nabla T[/math]. El gradiente de nuestro campo escalar indica la dirección en la que aumenta la temperatura, y su módulo indica cuánto aumenta.

Este gradiente se puede calcular mediante la siguiente fórmula: [math]\nabla T = \frac{\partial T}{\partial x}\vec i+\frac{\partial T}{\partial y}\vec j+\frac{\partial T}{\partial z}\vec k[/math] Con lo cual, obtenemos la siguiente expresión: [math]\nabla T = -2xsin(x^2)\vec i+2cos((y-1)^2)(y-1)\vec j[/math]

Para visualizar las curvas isotermas, se ha empleado el siguiente código.

Curvas de Nivel
Curvas de nivel de la función temperatura T

clear
clc
% Definir la placa
h=0.2;
r=1:h:2;
a=atan(1/2);
t=linspace(a,pi-a,pi/h);
[R,Th]=meshgrid(r,t);
% Conversión a cilíndricas
X=R.*cos(Th);
Y=R.*sin(Th);
% Campo de temparatura
T=cos(X.^2)+sin((Y-1).^2);
% Representación de las curvas de nivel
figure
mesh(X,Y,0*X);
hold on
contour(X,Y,T,20)
colorbar
title('CURVAS DE NIVEL DE LA TEMPERATURA')
axis ([-3, 3, -1, 3,])
view(0,90);
hold off
Tmax=max(max(T));

A partir de la imagen superior derecha, se puede ver que las temperaturas más altas se encuentran en la parte superior del sólido. En concreto, la temperatura máxima es de 1.8150 grados en el punto (0,2), obtenido en la línea 24 del código superior.

A continuación, habiendo calculado previamente el gradiente de la temperatura, se puede comprobar como el gradiente es ortogonal a las curvas de nivel.

Gradiente de la temperatura
Gradiente y curvas de nivel de la función temperatura T

clear
clc
% Definir la placa
h=0.2;
r=1:h:2;
a=atan(1/2);
t=linspace(a,pi-a,pi/h);
[R,Th]=meshgrid(r,t);
% Conversión a cilíndricas
X=R.*cos(Th);
Y=R.*sin(Th);
% Campo de temparatura
T=cos(X.^2)+sin((Y-1).^2);
% Representación de las curvas de nivel
figure
mesh(X,Y,0*X);
hold on
% Cálculo del gradiente
Gx=-2.*X.*sin(X.^2);
Gy=2.*cos((Y-1).^2).*(Y-1);
contour(X,Y,T,20)
quiver(X,Y,Gx,Gy);
colorbar
title('CURVAS DE NIVEL Y GRADIENTE')
axis ([-3, 3, -1, 3,])
view(0,90);
hold off

3 Campo vectorial de la energía calorífica

Según la Ley de Fourier, la energía calorífica [math]\vec Q[/math] viaja siguiendo la expresión: [math] \vec Q =-k▽T[/math], siendo k la constante de conductividad térmica del material. Para este ejemplo supondremos la constante de conductividad térmica de la placa como k=1, por lo que nuestra ecuación queda como: [math]\vec Q=-▽T[/math].

Tomando el gradiente del campo de temperaturas en coordenadas cartesianas: [math]▽T(x,y)=-2xsen(x^2)\vec i +2(y-1)cos((y-1)^2)\vec j [/math]:

[math]\vec Q =2xsen(x^2)\vec i -2(y-1)cos((y-1)^2)\vec j [/math]

Energía calorífica de la placa

clc
% Definición de variables
paso_muestreo = 2/10;
rho = 1:paso_muestreo:2;
ang = atan(1/2);
theta = linspace(ang,pi-ang,pi/paso_muestreo);
[D,A] = meshgrid(rho,theta);

% Conversión a cartesianas
X = D.*cos(A);
Y = D.*sin(A);
mesh(X,Y,0*X);

% Aplicación del campo vectorial a los puntos de la placa
Qx = 2.*X.*sin(X.^2);
Qy = 2.*(Y-1).*cos((Y-1).^2);
hold on
quiver(X,Y,Qx,Qy);
hold off

view(2)
axis([-3,3,-1,3])
% Definición del título
title('ENERGÍA CALORÍFICA DE LA PLACA')

4 Representación del campo de vectores en los puntos del sólido

Representamos el campo de vectores de deslizamiento en el cuarto de anillo, para eso nos ofrecen un campo que va a ser: [math] \vec u = \frac{1}{2}sen(θ) \vec e_ρ [/math]. Pasamos nuestro \vec u a cartesianas, sabiendo que [math] \vec e_ρ= cos(θ) \vec i + sin(θ)\vec j[/math]

clc
paso_muestreo = 2/10;
  distancia = 1:paso_muestreo:2;
  a = atan(1/2);
  angulo = linspace(a,pi-a,pi/paso_muestreo);
  [D,A] = meshgrid(distancia,angulo);
  % Conversión a cilíndricas
  x = D.*cos(A);
 y = D.*sin(A);
 mesh(x,y,0.*x)
 view(2)
 axis([-3,3,-1,3])
 hold on
 mx=0.5.*sin(A).*cos(A);
 my=0.5.*sin(A).*sin(A);
 quiver(x,y,mx,my);
 % Definición del título
 title('CAMPO DE VECTORES')

5 Imagen del sólido antes y después del desplazamiento

A continuación, se representará el sólido antes y después de la deformación producida por el campo de presiones [math] \vec u [/math] en coordenadas cartesianas (calculado en el apartado anterior). Se puede observar además que la deformación solo irá en la dirección [math] \vec e_{\rho} [/math] al ser [math] \vec e_{\theta} [/math] nula.

right

clc
% PLACA SIN DEFORMAR
paso_muestreo = 2/10;                         
rho = 1:paso_muestreo:2;
alfa = atan(1/2);
theta = linspace(alfa,pi-alfa,pi/paso_muestreo);
[RHO,THETA] = meshgrid(rho,theta);
x = RHO.*cos(THETA);
y = RHO.*sin(THETA);
subplot(1,2,1)                  
Q = mesh(x,y,0*x);
view(2)
set(Q,'EdgeColor','b');
axis ([-3,3,-1,3])
title('PLACA SIN DEFORMAR');
% PLACA DEFORMADA
subplot(1,2,2)                     
A = (sin(THETA).*cos(THETA))./2;
B = (sin(THETA).^2)./2;
X = x+A;
Y = y+B;
S = mesh(X,Y,0*X);
view(2)
set(S,'EdgeColor','r');
axis ([-3,3,-1,3])
title('PLACA DEFORMADA');
% COMPARACIÓN AMBAS PLACAS
hold on
figure                                 
S = mesh(X,Y,0*X);
axis ([-3,3,-1,3])
view(2)
set(S,'EdgeColor','r');
title('COMPARACIÓN DE AMBAS PLACAS');
hold on
Q = mesh(x,y,0*x);
set(Q,'EdgeColor','b');
hold off

center

6 Divergencia y sus mínimos y máximos

Para calcular la divergencia, empleamos la siguiente expresión: [math]\nabla ·u = \frac{1}{\rho} (\frac{\partial }{\partial \rho} (\rho u_\rho)+\frac{\partial }{\partial θ} (u_θ)+\frac{\partial }{\partial z} (\rho u_z))[/math]

Una vez operado, obtenemos

[math]\nabla ·u = \frac{senθ}{2\rho}[/math]

Con ello, representamos gráficamente la divergencia en Matlab

Divergencia
Representa la medida de cómo un campo se desplaza. Cuanto mayor el número, más se desplaza

clear
clc
% Definir la placa
h=0.2;
r=1:h:2;
a=atan(1/2);
t=linspace(a,pi-a,pi/h);
[R,Th]=meshgrid(r,t);
% Conversión a cartesianas
X=R.*cos(Th);
Y=R.*sin(Th);
% Cálculo de la divergencia
D=(sin(Th)./2./R);
surf(X,Y,D)
view(0,90);
axis equal
axis([-3,3,-1,3])
colorbar
title('DIVERGENCIA')
% Obtención de máximos y mínimos
Dmax=max(max(D));
Dmin=min(min(D));

A partir de la imagen superior derecha, se puede ver que las divergencias más altas se encuentran en la parte central del sólido. En concreto, la divergencia máxima es de 0,5 en el punto (0,1), y la mínima es de 0,1118 en los dos extremos longitudinales, obtenido en la línea 21 y 22 del código superior. Como la divergencia mínima es mayor a 0, no existe ningún punto de divergencia nula.

7 Rotacional en t = 0

Sea [math] ∇ × \vec u [/math] el rotacional de un campo de desplazamientos [math] \vec u = (ux, uy, uz)[/math] expresado en un sistema de coordenadas de vectores unitarios ortogonales [math]\vec{a}, \vec{b}, \vec{c}[/math], aplicado a una superficie bidimensional expresado en dicho sistema de coordenadas, se cumple que:

[math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{u} & \vec{v} &\vec{w} \\ d/da & d/db & d/dc\\ ux & uy & uz \end{vmatrix}[/math]

Continuando con nuestro ejemplo, calcularemos el rotacional del campo vectorial [math] \vec u (ρ,θ,z) = \frac{1}{2}sen(θ)·{\vec e_ρ} [/math] como:

[math] \mathbb{\nabla \times \vec u} = {1 \over ρ} \; \begin{vmatrix} {\vec e_ρ} & {ρ· \vec e_θ} & {\vec e_z}\\ {\partial \over \partial ρ} & {\partial \over \partial θ} & {\partial \over \partial z}\\ {\frac{1}{2}sen(θ)} & {0} & {0}\\ \end{vmatrix} [/math] = [math] -{1 \over ρ}·{\partial \over \partial θ}[{\frac{1}{2}sen(θ)}]·{\vec e_z} = [\frac{-1}{2} \frac{cos(θ)}{ρ}]·{\vec e_z}[/math]

Por lo tanto, [math] |∇ × \vec u |= \frac{1}{2} \frac{cos(θ)}{ρ} [/math]

Rotacional del campo vectorial en la placa

clc
% Definición de variables
paso_muestreo = 2/10;
rho = 1:paso_muestreo:2;
ang = atan(1/2);
theta = linspace(ang,pi-ang,pi/paso_muestreo);
[D,A] = meshgrid(rho,theta);

% Conversión a cartesianas
X = D.*cos(A);
Y = D.*sin(A);
mesh(X,Y,0*X);

% Definición del rotacional del campo
ROTu = 1/2.*cos(A)./D;

% Representación gráfica del rotacional
surf(X,Y,ROTu);
view(2)
axis([-3,3,-1,3])
colorbar
% Definición del título
title('REPRESENTACIÓN DEL ROTACIONAL DEL CAMPO')

Como podemos observar en la gráfica, los puntos con un mayor rotacinal son aquellos con menor ángulo θ, que darán un mayor valor al coseno; y con menor ρ.

8 Tensiones normales

Sea [math]\sigma=λ\bigtriangledown \cdot \vec{u}Ι + 2µε[/math] el campo de tensiones en un medio elástico, isótropo y homogéneo donde [math]ε=\frac{\bigtriangledown\vec{u}+\bigtriangledown\vec{u}^t}{2}[/math] la parte simétrica del gradiente del campo vectorial [math]\vec u[/math], los coeficientes de Lamé [math]λ,µ =1[/math] y [math]\bigtriangledown \cdot \vec{u}[/math] la divergencia del campo vectorial [math]\vec u[/math] se puede proceder a los cálculos de las deformacciones en las direcciones de los vectores de la base ortonormal orientada positivamente [math]\vec e_ρ,\vec e_θ,\vec e_z[/math]. Para ello, se hacen unos cálculos intermedios que facilitan la operación:

Calculamos el gradiente del campo y su traspuesto:

[math]\bigtriangledown\vec{u}=\begin{pmatrix} \ 0 & \frac{cos(θ)}{2\rho} & 0 \\ 0 & \frac{sen(θ)}{2\rho} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]\bigtriangledown\vec{u^t}=\begin{pmatrix} \ 0 & 0 & 0 \\ \frac{cos(θ)}{2\rho} & \frac{sen(θ)}{2\rho} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]ε=\begin{pmatrix} \ 0 & \frac{cos(θ)}{4\rho} & 0 \\ \frac{cos(θ)}{4\rho} & \frac{sen(θ)}{2\rho} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math] ––––––––––>[math]2ε=\begin{pmatrix} \ 0 & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{sen(θ)}{\rho} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

A continuación calculamos la divergencia que resulta igual a: [math]\nabla\cdot\vec u=\frac{sen(θ)}{2\rho}.[/math]

Por tanto ya podemos calcular σ:

[math] σ= \begin{pmatrix} \ \frac{sen(θ)}{2\rho} & 0 & 0 \\ 0 & \frac{sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix} + \begin{pmatrix} \ 0 & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{sen(θ)}{\rho} & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}[/math]

Por tanto, representaremos las tensiones normales en la dirección que marca el eje [math]\overrightarrow e_{\rho} [/math], es decir [math]\overrightarrow e_{\rho} \cdot σ \cdot \overrightarrow e_{\rho}[/math], las tensiones normales en la dirección que marca el eje[math]\overrightarrow e_{\theta} [/math], es decir [math]\overrightarrow e_{\theta} \cdot σ \cdot \overrightarrow e_{\theta}[/math] y las correspondientes al eje [math]\overrightarrow e_{z} [/math], es decir [math]\overrightarrow e_{z} \cdot σ \cdot \overrightarrow e_{z}[/math]. A continuación las calculamos y representamos:

En primer lugar, dirección de [math]\vec e_\rho[/math]-->

[math]\overrightarrow e_{\rho} \cdot σ \cdot \overrightarrow e_{\rho}=\begin{pmatrix} 1\\0\\0 \end{pmatrix}\cdot\begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}\cdot\begin{pmatrix} 1\\0\\0 \end{pmatrix}=\frac{sen(θ)}{2\rho}[/math]

En segundo lugar, dirección de [math]\vec e_\theta[/math]-->

[math]\overrightarrow e_{\theta} \cdot σ \cdot \overrightarrow e_{\theta}=\begin{pmatrix} 0\\1\\0 \end{pmatrix}\cdot\begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}\cdot\begin{pmatrix} 0\\1\\0 \end{pmatrix}=\frac{3sen(θ)}{2\rho}[/math]

Por último, dirección de [math]\vec e_z[/math]-->

[math]\overrightarrow e_{z} \cdot σ \cdot \overrightarrow e_{z}=\begin{pmatrix} 0\\0\\1 \end{pmatrix}\cdot\begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}\cdot\begin{pmatrix} 0\\0\\1 \end{pmatrix}=\frac{sen(θ)}{2\rho}[/math]

Para finalizar el apartado de tensiones, representaremos cada una de las direcciones que hemos calculado:

clc
% Definicion de variables
paso_muestreo = 2/10;                         
rho = 1:paso_muestreo:2;
alfa = atan(1/2);
theta = linspace(alfa,pi-alfa,pi/paso_muestreo);
[RHO,THETA] = meshgrid(rho,theta);
x = RHO.*cos(THETA);
y = RHO.*sin(THETA);

t1= 1/2.*sin(THETA)./RHO;
subplot(3,1,1)
surf(x,y,t1)
axis equal
title('Tensión normal en la dirección del eje rho')

t2= 3/2.*sin(THETA)./RHO;
subplot(3,1,2)
surf(x,y,t2)
axis equal
title('Tensión normal en la dirección del eje theta')

t3= 1/2.*sin(THETA)./RHO;
subplot(3,1,3)
surf(x,y,t3)
axis equal
title('Tensión normal en la dirección del eje z')

9 Tensiones tangenciales

Se calculará las tensiones tangenciales respecto al plano ortogonal a [math]\vec e_\rho[/math]. Se pide además las no nulas.

Del apartado anterior sacamos que [math] σ= \begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}[/math] y [math]\overrightarrow e_{\rho} \cdot σ \cdot \overrightarrow e_{\rho}=\frac{sen(θ)}{2\rho}[/math]. Por tanto:

[math]|σ·\vec e_ρ-(\vec e_ρ·σ·\vec e_ρ)·\vec e_ρ| = \left | \begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}·\begin{pmatrix} 1\\0\\0 \end{pmatrix} - \frac{sen(θ)}{2\rho}·\begin{pmatrix} 1\\0\\0 \end{pmatrix}\right | = \left | \begin{pmatrix} 0\\ \frac{cos(θ)}{2\rho}\\0 \end{pmatrix}\right| [/math]

Visualización 3D de las tensiones tangenciales

clc
% Definición de variables
paso_muestreo = 2/10;
rho = 1:paso_muestreo:2;
alfa = atan(1/2);
theta = linspace(alfa,pi-alfa,pi/paso_muestreo);
[RHO,THETA] = meshgrid(rho,theta);
x = RHO.*cos(THETA);
y = RHO.*sin(THETA);
% Tensiones tangenciales respecto al eje RHO
Sigma = cos(THETA)./(2.*RHO);
surf(x,y,Sigma);
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title('TENSIONES TANGENCIALES');

10 Tensión de Von Mises

Richard Von Mises

La tensión de Von Mises debe su nombre al ingeniero austrohúngaro Richard Von Mises, del cual podemos encontrar una foto a la derecha.

Es una magnitud que se utiliza en la mecánica de materiales y el análisis de estructuras, que se emplea para anticipar el fallo de un material que está soportando una carga, y que puede llegar a sufrir deformaciones permanentes. Por lo tanto, si la tensión de Von Mises supera el límite de fluencia, dejará de deformarse elásticamente para empezar a deformarse plásticamente.

Esta viene determinada por la siguiente fórmula en la que σ1, σ2 y σ3 son los autovalores del tensor de tensiones σ(ρ,θ).

[math] \sigma_V =\sqrt{\frac{(\sigma_1-\sigma_2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2}{2}}[/math]

Estos autovalores los obtenemos de la matriz obtenida en el apartado 8: [math] σ= \begin{pmatrix} \ \frac{sen(θ)}{2\rho} & \frac{cos(θ)}{2\rho} & 0 \\ \frac{cos(θ)}{2\rho} & \frac{3sen(θ)}{2\rho} & 0 \\ 0 & 0 & \frac{sen(θ)}{2\rho} \end{pmatrix}[/math]

Los autovalores calculados son: [math] σ_1=\frac{sen(θ)}{2\rho}[/math] , [math]σ_2=\frac{2sen(θ)-1}{2\rho}[/math] y [math]σ_3=\frac{2sen(θ)+1}{2\rho}[/math]

Sustituyendo estos valores en la fórmula de la tensión de Von Mises, obtenemos la siguiente expresión:

[math] \sigma_V =\sqrt{\frac{(\frac{sen(θ)}{2\rho}-\frac{2sen(θ)-1}{2\rho})^2+(\frac{2sen(θ)-1}{2\rho}-\frac{2sen(θ)+1}{2\rho})^2+(\frac{2sen(θ)+1}{2\rho}-\frac{sen(θ)}{2\rho})^2}{2}}=\frac{\sqrt{sen^2(θ)+3}}{2\rho}[/math]

Tensión de Von Mises
Límite de la deformación elástica

Para calcular las tensiones de Von Mises se ha utilizado el siguiente código, dando como resultado la imagen que se encuentra a la derecha. Además, la tensión máxima del anillo, calculada en la última línea del código es de 1MPa en el punto (0,1).

clear
clc
% Definir la placa
h=0.2;
r=1:h:2;
a=atan(1/2);
t=linspace(a,pi-a,pi/h);
[R,Th]=meshgrid(r,t);
% Conversión a cilíndricas
X=R.*cos(Th);
Y=R.*sin(Th);
% Tensión de Von Mises
VM=zeros(15,6);
VonMises=@ (t1,t2,t3) sqrt((((t1-t2)^2+(t2-t3)^2+(t3-t1)^2)/2));
% Cálculo de la tensión de Von Mises para cada punto del mallado
for i=1:15
    for j=1:6
        MTensiones=[sin(Th(i,j))./(2.*R(i,j)),cos(Th(i,j))./(2.*R(i,j)),0;cos(Th(i,j))./(2.*R(i,j)),3*sin(Th(i,j))./(2.*R(i,j)),0;0,0,sin(Th(i,j))./(2.*R(i,j))];
        Autov=eig(MTensiones);
        aut1=Autov(1);
        aut2=Autov(2);
        aut3=Autov(3);
        VM(i,j)=VonMises(aut1,aut2,aut3);
    end
end
surf(X,Y,VM)
title('TENSIÓN DE VON MISES')
axis ([-3, 3, -1, 3,])
view(0,90);
colorbar
mvm=max(max(VM));