Estudio de Transformaciones de un Cuarto de Anillo (43-C)

De MateWiki
Saltar a: navegación, buscar
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:

[math] \vec u(ρ,θ) = \frac{log(3 - ρ)}{2}cos(2θ)\vec e_ρ [/math]


Temperatura, T:

[math]T(x, y)=sin((x^2 + (y - 3)^2)[/math]


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.

right
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)

Gráfica Pregunta 2A.jpg

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

Por lo tanto el gradiente será:

[math]\nabla T(x,y) =cos(x^2+(y-3)^2)2x\vec i + cos(x^2+(y-3)^2)2(y-3)\vec j [/math]
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

GradienteDef.jpg 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:

Q=-k∇T


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 off

Fourier.jpg

4 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:

[math] \vec u(ρ,θ) = \frac{log(3-ρ)}{2}cos(2θ)\vec e_ρ [/math]


right
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 off
center
center

6 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:

derecha
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;

[math]\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec e_ρ & ρ\vec e_θ & \vec e_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right| [/math]

aplicando el campo de desplazamiento, devuelve la expresión;

[math]\nabla×\vec u(ρ,θ)=\frac{ln(3-p)(sen(2θ))}{ρ} \vec e_z [/math]

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.

derecha
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:

[math]\sigma = λ \nabla · \vec u I + 2µ \epsilon [/math]



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

Representación del mallado
Representación del mallado
%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:

[math] T = (\vec σ·\vec n - \vec n\otimes\vec n [/math]·σ)[math]\vec t[/math]



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]

right
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.

right
%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');