Estudio de Transformaciones de un Cuarto de Anillo (43-C)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Estudio de Transformaciones de un Cuarto de Anillo (43-C) |
| Asignatura | Teoría de Campos |
| Curso | 2023-24 |
| Autores | Samantha Ugarte Flores, Celia Bolívar Illana, Miguel Goñi Fernández de Aguirre, Enrique Echevarría del Alcázar |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Se considera una placa plana que ocupa un cuarto de un anillo circular centrado en el origen y comprendido entre los radios 1 y 2, en el plano y ≥|x|/2. Dibujar con los ejes en el cuadrado [-3,3] x [-1,3]. Trabajando en coordenadas cilíndricas.
Con el siguiente campo:
Temperatura, T:
Contenido
1 Mallado
El sólido está definido por dos condiciones. Por un lado, el radio está contenido entre 1 y 2. Además el sólido está contenido en el plano y ≥|x|/2. Esta última condición se traduce en coordenadas cilíndricas en que rsenθ≥|rcosθ/2|. Simplificando esta expresión se deduce que tanθ≥1/2 y tanθ<-1/2. Los ángulos que cumplen está condición son los comprendidos entre (atan(1/2),pi-atan(1/2)) por lo que θ se mueve entre dichos valores. El código que aparece a continuación es el que permite representar en Matlab el sólido definido por estas condiciones.
h=0.1; %paso de muestreo
r=1:h:2; %condición de ρ
tt=linspace(atan(1/2),pi-atan(1/2),pi*10); %condición de θ
[RR,TT]=meshgrid(r,tt); %mallado
x=RR.*cos(TT);
y=RR.*sin(TT);
mesh(x,y,0.*x)
view(2)
axis ([-3,3,-1,3]) %restricción de ejes
xlabel('Eje X')
ylabel('Eje Y')
title('Representación en 2D de la placa plana')
2 Curvas de nivel y gradiente
La temperatura viene dada por la siguiente función:
[math]T(x, y)=sin((x^2 + (y - 3)^2)[/math]
En el siguiente gráfico está representado a la izquierda las curvas de nivel de la función de temperatura y a la derecha hay una representación en 3D de la función.
h= 0.2; %Paso de muestreo
r= 1:h:2;
t= linspace(atan(1/2),pi-atan(1/2),pi*10);
[rr,tt]= meshgrid(r,t); %Mallado
xx=rr.*cos(tt);
yy=rr.*sin(tt);
T=sin(xx.^2+(yy-3).^2);
hold on
subplot(1,2,2) %Dividimos la pantalla en dos
surf(xx,yy,T); %representación de la temperatura en 3D
colorbar;
axis([-3,3,-1,3]);
axis vis3d
title('Temperatura en 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje z')
subplot(1,2,1) %Escribimos en la segunda pantalla
contour(xx,yy,T,40) %Líneas de nivel
title('Curvas de Nivel','Fontsize',16)
colorbar %Mostramos las escala
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
hold off
maximo=max(max(T))
fprintf('La máxima temperatura es: %1.4f \n',T)Por lo tanto el gradiente será:
h =0.2; %Paso de muestreo
r = 1:h:2;
t=linspace(atan(1/2),pi-atan(1/2),35);
[rr,tt]= meshgrid(r,t); %Mallado y parametrización
xx=rr.*cos(tt);
yy=rr.*sin(tt);
Gradx=cos(xx.^2+(yy-3).^2.*2.*xx) ;
Grady=cos(xx.^2+(yy-3).^2.*2.*(yy-3)) ; %GRADIENTE
Temperatura=sin(xx.^2+(yy-3).^2); %TEMPERATURAS
subplot(1,2,2) %Creación del gradiente
view(2)
hold on
surf(xx,yy,Temperatura);
axis([-3,3,-1,3]);
quiver3(xx,yy,Temperatura,Gradx,Grady,0*Temperatura); %Dibujo del gradiente
axis([-3,3,-1,3]);
axis vis3d
view(3)
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Gradiente 3D')
subplot(1,2,1)
mesh(xx,yy,0.*xx)
view(2)
axis([-3,3,-3,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Gradiente 2D')
hold on
quiver(xx,yy,Gradx,Grady)
hold off
Como se puede observar, el gradiente es perpendicular a la superficie.
3 Ley de Fourier
De acuerdo con la Ley de Fourier la energía calorífica Q es:
donde k es la constante de conductividad térmica de la placa que supondremos k=1.
Por ello la conductividad térmica será el gradiente de la temperatura (calculado en el apartado anterior) pero en sentido contrario (en negativo).
h=2/10;
r=1:h:2;
t=pi/8:h:7*pi/8;
[U,T]=meshgrid(r,t);
x=U.*cos(T);
y=U.*sin(T);
Temp=sin(x.^2+(y-3).^2);
Tempx=-(2*x.*cos(x.^2+(y-3).^2));
Tempy=-((2*y-6).*cos(x.^2+(y-3).^2));
contour(x,y,Temp);
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar;
hold on
quiver(x,y,Tempx,Tempy);
axis equal
axis([-3,3,-1,3]);
view(2);
title('Gradiente del campo T');
hold off4 Campo de vectores
En este apartado mostraremos el campo de vectores en los puntos del mallado del solido. Nos piden los resultados en t=0, pero nuestro campo no depende de t.
Nuestro campo es:
h=2/10;
r=1:h:2;
t=pi/4:h:3.*pi/4;
[U,T]=meshgrid(r,t);
x=U.*cos(T);
y=U.*sin(T);
FX=((log(3-U)/2).*cos(2*T).*cos(T));
FY=((log(3-U)/2).*cos(2*T).*sin(T));
quiver(x,y,FX,FY);
axis ([-3,3,-1,3]);
xlabel('Eje X')
ylabel('Eje Y')
view(2)
axis equal
title('Campo de vectores')
5 Desplazamiento de la placa
Estudiamos el desplazamiento que sufre el cuarto de anillo gracias al campo [math] \vec u(ρ,θ) [/math].
A continuación hemos puesto los gráficos que representan este desplazamiento y su código.
h=2/10; %Paso de muestreo
rr=1:h:2; %Asignación de las variables rho y theta
tt=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,TT]=meshgrid(rr,tt); %Mallado
x=RR.*cos(TT); %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
subplot(1,2,1) %Placa en el inicio
m=mesh(x,y,0*x);
view(2)
set(m,'EdgeColor','b') %Color de la placa de inicio
axis([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Placa no desplazada')
subplot(1,2,2) %Placa desplazada
X=((log(3-RR)/2).*cos(2.*TT)).*cos(TT)+x;
Y=((log(3-RR)/2).*cos(2.*TT)).*sin(TT)+y;
n=mesh(X,Y,0*X);
view(2)
set(n,'EdgeColor','r') %Color de la placa desplazada
axis([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Placa desplazada')
figure %Comparación de ambas placas
n=mesh(X,Y,0*X);
axis([-3,3,-1,3])
view(2)
set(n,'EdgeColor','r') %Color de la placa desplazada
title('Desplazamiento de la placa')
hold on
m=mesh(x,y,0*x);
set(m,'EdgeColor','b') %Color de la placa de inicio
hold off6 Divergencia
La divergencia de un campo vectorial mide el cambio de volumen debido al desplazamiento que sufre.
El campo que tenemos es [math] \vec u(ρ,θ) = \frac{log(3-ρ)}{2}cos(2θ)\vec e_ρ [/math] ,que con la fórmula de la divergencia en cilíndricas [math]\nabla\cdot\vec u=\frac{1}{ρ}[\frac{\partial}{\partial{ρ}}({ρ·u_ρ})] [/math] nos da:
[math]\nabla\cdot\vec u= \frac{-cos(2θ)}{2(3-ρ)}+\frac{log(3-ρ)cos(2θ)}{2ρ}[/math]
6.1 Gráfico
Teniendo la divergencia, hemos sacado su gráfico con el siguiente código:
h=2/10; %Paso de Muestreo
rr=1:h:2; %Asignación de las variables rr y tt
tt=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,TT]=meshgrid(rr,tt); %Mallado
x=RR.*cos(TT); %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
DIVu=(-cos(2.*TT))./(2.*(3-RR))+(log(3-RR).*cos(2.*TT)./2.*RR); %Divergencia de u
subplot(1,2,1) %Gráfica 2D
surf(x,y,DIVu)
shading faceted
view(2)
axis([-3,3,-1,3])
colorbar
title('Divergencia 2D')
xlabel('Eje X')
ylabel('Eje Y')
subplot(1,2,2) %Gráfica 3D
view(2)
surf(x,y,DIVu)
shading faceted
axis([-3,3,-1,3])
colorbar
axis vis3d
title('Divergencia 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
shg
6.2 Máx, mín, nulo
El punto con máxima divergencia es (0;2) con 0.5000 de valor; los puntos con mínima divergencia son (1.78885;0.894427) y (-1.78885;0.894427) con -0.3000. La divergencia es nula cuando θ=π/4.
DIVumax=max(max(DIVu));
DIVumin=min(min(DIVu));
7 Rotacional
El rotacional es un operador vectorial que, muestra la tendencia de un campo vectorial a producir rotación en un punto. La expresión del rotacional en coordenadas cilíndricas viene definida por la siguiente expresión;
aplicando el campo de desplazamiento, devuelve la expresión;
En las gráficas se puede observar como el punto que sufre un mayor rotacional es el (0,1).Dichas gráficas han sido creadas con el siguiente código.
h=0.1; %Paso de muestreo
rr=1:h:2; %Usamos coordenadas polares
tt=linspace(atan(1/2),pi-atan(1/2),16);
[RR,TT]= meshgrid(rr,tt); %Mallado
xx=RR.*cos(TT);
yy=RR.*sin(TT);
ROT=log(3-RR)./RR.*sin(2.*TT); %Modulo del rotacional de u
subplot(1,2,1); %Rotacional en 2D
surf (xx,yy,ROT);
axis([-3,3,-1,3]);
view(2);
colorbar;
title('Rotacional en 2D');
xlabel('Eje X')
ylabel('Eje Y') %Rotacional en 3D
subplot(1,2,2);
surf(xx,yy,ROT);
colorbar
title('Rotacional en 3D');
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
8 Tensiones
Los desplazamientos en un medio elástico lineal, isótropo y homogéneo describen el tensor de tensiones con la fórmula:
Siendo [math]\epsilon [/math] la parte simétrica del tensor gradiente de [math]\vec u [/math] conocido como tensor de deformaciones.
Mientras que λ , μ son los coeficiente de Lamé, los cuales dependen de las cualidades elásticas del material. Los dos parámetros juntos constituyen una parametrización del módulo de elasticidad para medios isótropos homogéneos. En este caso las propiedades de la placa son tales que resultan ambos coeficientes ser 1.
Después podremos representar las tensiones normales en la dirección que marca el eje i, es decir [math]\vec i [/math] = \(\vec g_{i}\) * σ * \(\vec g_{i}\), las tensiones normales en la dirección que marca el eje j, es decir [math]\vec j [/math] = \(\vec g_{j}\) * σ * \(\vec g_{j}\) y las correspondientes al eje k, es decir [math]\vec k[/math] = \(\vec g_{k}\) * σ * \(\vec g_{k}\).
8.1 Tensor de Deformaciones
Se define el tensor de deformaciones [math]\epsilon(\vec u(ρ,θ))[/math] como la parte simétrica del tensor [math] \nabla{\vec u(ρ,θ)} [/math],
- [math]\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}[/math].
Para ello se calculan [math]\nabla{\vec u(ρ,θ)}[/math] y [math](\nabla{\vec u(ρ,θ)})^t[/math].
- [math]\frac{\partial \vec u}{\partial ρ}=(\frac{-cos(2θ)}{6-2ρ}) \vec e_p[/math]
- [math]\frac{\partial \vec u}{\partial θ}= -log(3-ρ)sen(2θ) \vec e_p + (\frac{ln(3-ρ)}{2})cos(2θ) \vec e_θ[/math]
- [math]\frac{\partial \vec u}{\partial z}= 0 \vec e_z[/math]
Por lo tanto el tensor de deformaciones será igual a:
- [math]\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} (\frac{-cos(2θ)}{6-2p}) & (\frac{-ln(3-p)(sen(2θ))}{2}) & 0 \\ (\frac{-ln(3-p)(sen(2θ))}{2}) & (ln(3-ρ)cos(2θ) & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].
8.2 Tensor de Tensiones
El tensor de tensiones se define con la divergencia, el tensor de deformaciones y los coeficientes de Lamé. Por tanto, siguiendo su expresión queda:
[math]\sigma = \left(\begin{matrix} (\frac{ln(3-p)cos(2θ)}{2ρ})-\frac{cos2θ)}{6-2p}) & 0 & 0 \\ & \\0 & (\frac{ln(3-p)cos(2θ)}{2ρ})-\frac{cos2θ)}{6-2p}) & 0 \\ & \\0 & 0 & (\frac{ln(3-p)cos(2θ)}{2ρ})-\frac{cos2θ)}{6-2p}) \end{matrix} \right) + 2 ·\left ( \begin{matrix} (\frac{-cos(2θ)}{6-2p}) & (\frac{-ln(3-p)(sen(2θ))}{2}) & 0 \\ (\frac{-ln(3-p)(sen(2θ))}{2}) & (ln(3-ρ)cos(2θ) & 0 \\ 0 & 0 & 0 \end{matrix} \right) [/math]
[math] = \left ( \begin{matrix} (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) & (-ln(3-p)(sen(2θ))) & 0 \\ (-ln(3-p)(sen(2θ))) & (\frac{(4p+1)ln(3-p)cos(2θ)}{2ρ}) & 0 \\ & \\0 & 0 & (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{cos(2θ)}{6-2p}) \end{matrix} \right )
[/math]
8.3 Tensiones Normales
A pesar de que el desplazamiento sea plano ([math]\vec u(θ) [/math] solo dependa de una dirrección) las tensiones pueden no serlo y pueden haber tensiones ortogonales al plano de la placa. Utilizando el siguiente tensor podremos comprobarlo:
-
[math]T = \vec n\otimes\vec n · \vec σ[/math]
-
8.3.1 Tensión normal en la dirección que marca el eje [math]\vec e_ρ [/math]
El módulo de la tensión normal es el siguiente:
[math] |(σ·\vec e_ρ))\vec e_ρ | = (\left ( \begin{matrix} (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) & (-ln(3-p)(sen(2θ))) & 0 \\ (-ln(3-p)(sen(2θ))) & (\frac{(4p+1)ln(3-p)cos(2θ)}{2ρ}) & 0 \\ & \\0 & 0 & (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{cos(2θ)}{6-2p}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right )) ·\vec e_ρ[/math]
[math] = (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) [/math]
8.3.2 Tensión normal en la dirección que marca el eje [math]\vec e_θ [/math]
El módulo de la tensión normal es el siguiente:
[math]|(σ·\vec e_θ))\vec e_θ | = (\left ( \begin{matrix} (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) & (-ln(3-p)(sen(2θ))) & 0 \\ (-ln(3-p)(sen(2θ))) & (\frac{(4p+1)ln(3-p)cos(2θ)}{2ρ}) & 0 \\ & \\0 & 0 & (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{cos(2θ)}{6-2p}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right )) ·\vec e_θ[/math]
[math] = (\frac{(4p+1)ln(3-p)cos(2θ)}{2ρ}) [/math]
8.3.3 Tensión normal en la dirección que marca el eje [math]\vec e_z [/math]
El módulo de la tensión normal es el siguiente:
[math]|(σ·\vec e_θ))\vec e_θ | = (\left ( \begin{matrix} (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) & (-ln(3-p)(sen(2θ))) & 0 \\ (-ln(3-p)(sen(2θ))) & (\frac{(4p+1)ln(3-p)cos(2θ)}{2ρ}) & 0 \\ & \\0 & 0 & (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{cos(2θ)}{6-2p}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\0\\1\end{matrix} \right )) ·\vec e_θ[/math]
[math] = (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{cos(2θ)}{6-2p}) [/math]
A partir de el siguiente código podemos observar las gráficas de los tensores
%Realizamos el mallado
r=1:0.2:2;
t=linspace(atan(1/2),pi-atan(1/2),10);
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
%Elemento (1,1,1) de la matriz sigma
M11=((log(3-rr).*(cos(2.*tt)))./(2.*rr)-(3.*cos(2.*tt)./(6-2.*rr)));
%Elemento (2,2,2) de la matriz sigma
M22=((4.*rr+1).*log(3-rr).*cos(2.*tt))./(2.*rr);
%Elemento (2,2,2) de la matriz sigma
M33=((log(3-rr).*(cos(2.*tt)))./(2.*rr)-(cos(2.*tt)./(6-2.*rr)));
subplot(1,3,1)
surf(xx,yy,M11)
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(2)
colorbar
title('Tensión normal en la dirección del eje e sub rho');
subplot(1,3,2)
surf(xx,yy,M22)
axis equal
view(2)
colorbar
xlabel('x')
ylabel('y')
zlabel('z')
title('Tensión normal en la dirección del eje e sub theta');
subplot(1,3,3)
surf(xx,yy,M33)
axis equal
view(2)
colorbar
xlabel('x')
ylabel('y')
zlabel('z')
title('Tensión normal en la dirección del eje e sub z');
figure
subplot(1,3,1)
surf(xx,yy,M11)
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
colorbar
title('Tensión normal en la dirección eje e sub rho');
subplot(1,3,2)
surf(xx,yy,M22)
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
colorbar
title('Tensión normal en la dirección eje e sub theta');
subplot(1,3,3)
surf(xx,yy,M33)
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
colorbar
title('Tensión normal en la dirección eje e sub z');
9 Tensiones tangenciales
Las tensiones en la dirección tangencial al plano de nuestra placa es la proyección de la tensión [math]\vec σ [/math] al plano, que tiene dirección: [math]\vec t[/math].
La proyección no es el tensor:
9.1 Tensión tangencial respecto al plano ortogonal al eje [math]\vec e_ρ [/math]
El módulo de la tensión tangencial es el siguiente:
- |[math]σ·\vec e_ρ-(\vec e_ρ·(σ·\vec e_ρ))\vec e_ρ | = [/math]
[math] \left | \left ( \begin{matrix} (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) & (-ln(3-p)(sen(2θ))) & 0 \\ (-ln(3-p)(sen(2θ))) & (\frac{(4p+1)ln(3-p)cos(2θ)}{2ρ}) & 0 \\ & \\0 & 0 & (\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{cos(2θ)}{6-2p}) \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right )-
(\frac{log(3-ρ)cos(2θ)}{2ρ})-\frac{3cos(2θ)}{6-2p}) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) \right | [/math]
[math] = (-ln(3-p)(sen(2θ))) [/math]
h=2/10; %Paso de muestreo
rho=1:h:2; %Asignación de las variables rho y theta
tt=linspace(atan(1/2),pi-atan(1/2),35);
[RR,TT]=meshgrid(rho,tt); %Mallado
x=RR.*cos(TT); %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
tang=((-(log(3-RR)).*(sin(2.*TT)))); %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 sub rho')
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 sub rho')
10 Ley de Von Mises
La tensión de Von Mises es una magnitud física escalar que indica el inicio del comportamiento plástico de un material. En ingeniería estructural se usa en el contexto de las teorías de fallo como indicador de un buen diseño para materiales dúctiles. La tensión de Von Misses se define por la fórmula:
[math]\sigma_{VM} = \sqrt {\frac{(\sigma_1 - \sigma_2 )^2 + (\sigma_2 - \sigma_3 )^2 + (\sigma_3 - \sigma_1 )^2}{2}}[/math], siendo [math] \sigma_1 , \sigma_2 , \sigma_3 [/math] autovalores de la tensión previamente calculada.
%Realizamos el mallado
r=1:0.2:2;
t=linspace(pi/4,3*pi/4,8);
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
%Calculamos la matriz de tensiones "sigma"
M11 = @(R,T) (log(R)./(R)).*cos(2.*(T)- (pi/2));
M12 = @(R,T) (sin(2.*T-(pi/2))./2.*(R)).*(1-log(R));
M22 = @(R,T) (3.*log(R)./(R)).*cos(2.*T-(pi/2));
M21 = @(R,T) (sin(2.*T-(pi/2))./(2.*(R))).*(1-log(R));
sigma = [];
VonMises = zeros(length(t), length(r));
for i = 1:length(t)
for j = 1:length(r)
sigma(1,1) = M11(rr(i,j),tt(i,j));
sigma(1,2) = M12(rr(i,j),tt(i,j));
sigma(1,3) = 0;
sigma(2,1) = M21(rr(i,j),tt(i,j));
sigma(2,2) = M22(rr(i,j),tt(i,j));
sigma(2,3) = 0;
sigma(3,1) = 0;
sigma(3,2) = 0;
sigma(3,3) = M11(rr(i,j),tt(i,j));
[f,c] = eig(sigma);
VonMises(i,j) = sqrt(((c(1,1)-c(2,2))^2+(c(2,2)-c(3,3))^2+(c(3,3)-c(1,1))^2)/2);
end
end
subplot(1,3,2);
surf(xx,yy,VonMises)
axis([-3,3,-1,3])
axis vis3d
colorbar;
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
title('Tensión de Von Mises');
subplot(1,3,3)
hold on
surf(xx,yy,VonMises);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');
MAX = max(VonMises,[],'all');
xlabel('Eje x')
zlabel('Eje z')
txt = ['Máxima tensión: ' num2str(MAX) ];
text(-2,0,0.87,txt)
hold off
subplot(1,3,1)
surf(xx,yy,VonMises);
view(2)
xlabel('Eje x')
ylabel('Eje y')
axis([-3,3,-1,3])
colorbar;
title('XOY');
