MateWiki - Contribuciones del usuario [es]

Visualización de campos escalares y vectoriales en elasticidad (Grupo 22-B)

2022-12-08T16:04:37Z

Sergio puente:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández}}
{{ Beta }}
== Introducción ==
El presente artículo tiene por objeto estudiar el desplazamiento experimentado por una pieza tras la aplicación de una fuerza.
 Para ello, se tienen dos cantidades físicas dependientes de las variables ρ y θ:
* La temperatura T(ρ,θ)
* El campo de desplazamientos <math>\vec u(ρ,θ)</math>, producido por la acción de una fuerza.
Para ello, se hará uso del software de programación y cálculo numérico Matlab/Octave.

== Descripción de la placa ==
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa de espesor diferencial con forma de cuarto de anillo circular, centrada en el origen y comprendida entre los radios 1 y 2 del plano <math>y≥|x|</math>, es decir, <math>1≤x^2+y^2≤4</math>.
Como se va a trabajar en coordenadas cilíndricas, esta región equivale a: <math>(ρ,θ)=[1,2]×[0,π/2]</math>. 
A continuación se muestra la sección de la placa.
[[Archivo:Placa plana.png|centre|]]
 
El código empleado para dibujar la placa es el siguiente:
{{matlab|codigo=

h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
mesh(x,y,z); %Dibujamos la placa
axis equal
axis([-3,3,-1,3]); %Ejes
view(2); %Vista en planta
title("Placa plana") %Nombre de la figura
}}

== Temperatura ==
La temperatura viene determinada por la siguiente función: <math> T(x,y)=log((x-3)^2)+2)</math>.
 En primer lugar se dibujarán las curvas de nivel de la temperatura.
[[Archivo:Curvas de nivel del campo de temperaturas.png|centre|]]
 El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel
axis equal;
axis([-3,3,-1,3]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
colorbar
}}
Ahora se dibujará la variación de la temperatura a lo largo de la placa.
[[Archivo:Campo de temperaturas.png|centre|]]
 '''El programa nos indica que la temperatura máxima alcanzada en la placa es de 3,0244ºC.'''
 El código empleado para dibujar el campo de temperatura a lo largo de la placa, así como para calcular la temperatura máxima alcanzada en la misma ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
surf(x,y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)
}}

== Gradiente de temperatura ==
Como se ha mencionado con anterioridad, el campo escalar que gobierna la temperatura en la placa es <math> T(x,y)=log((x-3)^2)+2)</math>.
 
El gradiente de un campo vectorial en coordenadas cartesianas se define como: 
<math>\nabla T(x,y,z) =\frac{d∂}{dx} +\frac{d∂}{dy} + \frac{d∂}{dz}</math>

De lo que resulta: <math>\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=>\frac{2(x-3)}{(x-3)^2+2} \vec i</math>
 
A continuación se muestra una representación del gradiente de la temperatura, como se puede apreciar este es ortogonal a las curvas de nivel de la temperatura ya que marca la dirección de máximo crecimiento de la misma.
Para este caso se observa que la temperatura crece en la dirección -<math>\vec i</math>
[[Archivo:Gradiente de temperatura_1.png|500px||center|]]
 
El código empleado para representar el gradiente ha sido el siguiente:
{{matlab|codigo=
h=0.2; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;;
T=log(((x-3).^2)+2); %Función de temperatura
grad=(2.*(x-3))./((x-3).^2+2); %Gradiente de temperatura
hold on
quiver(x,y,grad,0.*grad); %Se dibuja el gradiente
contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
axis([-3,3,-1,3]);
title("Gradiente de temperatura");
hold off
}}

== Campo de vectores ==
Se va a realizar una representación de el campo de desplazamientos que viene gobernado por la expresión: <math> \vec r(ρ,θ)= \frac{log(ρ)}{2}(sin(2θ-\frac{π}{2}))</math>
[[Archivo:Campo de desplazamientos sobre la curva.png|600px||center|]]
 
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=1/4*pi:pi/50:3/4*pi; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*(-sin(theta)); %Paso de la coordenada x del campo a cartesianas
FY=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*cos(theta); %Paso de la coordenada y del campo a cartesianas
quiver(X,Y,FX,FY); %Dibuja el campo vectorial
axis equal
axis([-3,3,-1,3]);
view(2) %Vista en planta
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')
}}
 

== Comparación de la placa antes y después del desplazamiento ==
En este apartado se va a realizar una comparación visual de la placa antes y después de sufrir el desplazamiento.
[[Archivo:Placa antes y después.png|800px||center|]]
[[Archivo:Comparación de placas.png|600px||center|]] 
El código empleado para llevar a cabo la comparación se expone a continuación:
{{matlab|codigo=
h=0.2;%paso de muestreo
r=1:h:2;
t=pi/4:h:((3*pi/4)+0.2);
[ro,theta]=meshgrid(r,t); %creamos el mallado
x=ro.*cos(theta);
y=ro.*sin(theta);
%sólido antes del desplazamieneto
subplot(1,2,1);
i=mesh(x,y,x.*0);
view(2);
set(i,'EdgeColor','r')
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa no desplazada');
%sólido despues del desplazamiento
subplot(1,2,2);
%pasamos el campo de desplazamiento a cartesianas
A=-(log10(ro)./2).*(sin(2.*theta-(pi/2)).*(sin(theta)));
B=(log10(ro)./2).*(cos(theta)).*(sin(2.*theta-(pi/2)));
%Este campo de desplazamiento no desplaza prácticamente la placa, si no que concentra la masa de la misma en el centroide.
X=x+A;
Y=y+B;
j=mesh(X,Y,0.*X);
view(2);
set(j,'EdgeColor','b');
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa desplazada');
%Comparación dde ambas placas
hold on
figure
j=mesh(X,Y,X*0);
axis equal
axis([-3,3,-1,3]);
view(2);
set(j,'EdgeColor','b');
title('Comparación de ambas placas: Roja (inicial), azul (final)');
hold on
i=mesh(x,y,x.*0);
set(i,'EdgeColor','r')
hold off
}}

==Divergencia==
La divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento. 
Para este supuesto se tiene que el campo de desplazamientos es <math> \vec u(ρ,θ) = \frac{ln(ρ)}{2}sin(2θ-\frac{π}{2})\vec e_θ </math> 
La divergencia resulta: <math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho·0})+\frac{\partial}{\partial{\theta}}(\frac{ln(ρ)}{2}sin(2θ-\frac{π}{2}))] </math>, es decir, <math>\nabla\cdot\vec u=\frac{ln(ρ)}{ρ}sin(2θ)</math>
 
La representación gráfica de la divergencia del campo de desplazamientos resulta: 
[[Archivo:divergencia_3.png|700px||center|]] 
El la placa, el valor máximo de la divergencia del campo de desplazamientos es 0.3466 y se alcanza en el punto (1.41421,1.41421). 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
%Paso de muestreo y definición de variables
h= 0.1;
r= 1:h:2;
t=pi/4:h:3*pi/4;
%Mallado
[rr,tt]= meshgrid(r,t);
%Parametrización
x=rr.*cos(tt);
y=rr.*sin(tt);
%Divergencia
div=(1./rr).*(log(rr).*sin(2.*(tt)));
%Divergencia en 2D
subplot(1,2,1)
surf(x,y,div)
shading interp
title('Divergencia en 2D')
axis([-3,3,-1,3]);
axis equal
colorbar;
view(2)
%Divergencia en 3D
subplot(1,2,2)
surf(x,y,div)
shading interp
axis([-3,3,-1,3]);
axis equal
colorbar;
view(3)
title('Divergencia en 3D')
maximo=max(max(div))
fprintf('El valor máximo de la divergencia es: %1.4f \n',maximo)
}}

== Rotacional ==
El rotacional de un campo vectorial permite observar la tendencia que tiene dicho campo a inducir rotación alrededor de un punto de la placa.
El rotacional, calculado en coordenadas cilíndricas viene determinado por la siguiente expresión:
 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|</math>. 
Si se aplica al campo de desplazamientos resulta: 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ·\frac{log(ρ}{2}·sin(2θ-\frac{π}{2}) & 0 \end{matrix}\right|=> -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}</math>. 
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
[[Archivo:Rotacional (2).png]] 
Como puede observarse, el rotacional es mayor en los puntos centrales de la placa, es por ello por lo que en la vista en 3 dimensiones la placa resulta doblada por esa zona. 
Del mismo modo, el programa nos indica que el valor máximo que alcanza el rotacional es 0.1504. 
El código empleado para dibujar el rotacional del campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
%Se crean los intervalos de las variables
h=50;
r=linspace(1,2,h);
t=linspace(pi/4,3*pi/4,h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=-(cos(2.*theta)./(2.*ro)).*log10((ro)+1)
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
shading interp
axis equal;
view(2);
colorbar;
title('Rotacional en 2D')
xlabel('Eje x');
ylabel('Eje y');
%Gráfica en 3D
subplot(1,2,2);
surf(x,y,rotacional);
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
Maximo=max(max(rotacional));
fprintf('el valor del rotacional máximo alcanzado es de: %1.4f \n',Maximo);
}}
==Tensor de deformaciones==
Las deformaciones que experimenta la placa se pueden modelizar mediante un tensor conocido como tensor de deformaciones <math>\epsilon(\vec u(ρ,θ))</math> que viene determinado por la parte simétrica del tensor <math> \nabla{\vec u(ρ,θ)} </math>, también llamado gradiente del campo de desplazamientos <math>\vec u(ρ,θ) </math>:
<math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> 
<math>\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}</math>
 
Lo primero será calcular <math>\nabla{\vec u(ρ,θ)}</math> que resulta: 
<math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial ρ} & \frac{\partial u_{1} }{\partial θ} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial ρ} & \frac{\partial u_{2} }{\partial θ} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial ρ} & \frac{\partial u_{3} }{\partial θ} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> 
 
 
::*<math>\frac{\partial \vec u}{\partial ρ}=(\frac{1}{2ρ})sen(2θ-\frac{π}{2}) \vec e_θ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial θ}= log(ρ)cos(2θ-\frac{π}{2}) \vec e_θ - (\frac{log(ρ)}{2})sen(2θ-\frac{π}{2}) \vec e_ρ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial z}= 0 \vec e_z</math>
 
 
El gradiente del campo de deformaciones resulta: 
<math>\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2})& 0\\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
A continuación se calcula la traspuesta de la matriz del gradiente del campo de deformaciones: 
<math>(\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
 
Aplicando la fórmula <math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> el tensor de deformaciones será igual a: 

<math>\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>

 
 

==Tensor de tensiones==
El tensor de tensiones viede deteminado por la siguiente expresión: <math>σ=λ·\nabla·\vec u I+2με</math> 
Desarrollando:
 
<math>\sigma = \left(\begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 & 0 \\ & \\0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) + 2 ·\left ( \begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix} \right) </math>
 
 
 
Finalmente resulta: <math>σ</math> = <math> \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right )
</math>
===Tensiones normales===
====Tensiones normales en la dirección del eje <math>\vec e_ρ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_ρ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math> 
 
Que resulta: <math>(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e rho.png|500px||left|]] [[Archivo:3drho_1.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
a=(log10(ro).*sin(2.*theta))./(ro);%elemento(1,1,1);
figure(1);
surf(x,y,a);
shading interp
axis equal
title('Tensión normal en la dirección e_ρ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_θ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_θ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math>
 
 
Que resulta: <math> (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e theta.png|500px|left]] [[Archivo:3drho.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
b=((log10(ro))./(ro))+2.*(log10(ro)).*sin(2.*theta);%elemento(2,2,2)
figure(1);
surf(x,y,b);
shading interp
axis equal
title('Tensión normal en la dirección e_θ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_z</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_z</math> es el siguiente: 
<math>|(σ·\vec e_θ))\vec e_θ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right ) ·\vec e_θ</math>
 
 
Que resulta: <math> (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
[[Archivo:Tension normal en e z.png|500px|left]] [[Archivo:3dz.png|500px|right]] 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha siso el siguiente: 
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
c=(log10(ro).*sin(2.*theta))./(ro);%elemento(3,3,3)
figure(1)
surf(x,y,c);
shading interp
axis equal
title('Tensión normal en la dirección e_k')
view(2);
colorbar;
}}

===Tensiones tangenciales respeto al plano ortogonal a <math>\vec e_ρ</math>===
Se van a calcular las tensiones tangenciales que sufre la placa en la dirección del eje <math>\vec e_ρ</math> en t=0. 
Para ello se conoce su el módulo es el siguiente: 
 
|<math>σ·\vec e_ρ-(\vec e_ρ·(σ·\vec e_ρ))\vec e_ρ ) </math>|. 
 
Desarrollando la expresión resulta: 
 
<math> \left | \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right )-
(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) \right | </math> 
 
Finalmente se obtiene: 
 
<math> \frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) </math>
La representación de la tensión tangencial que sufre la placa es la siguiente: 
[[Archivo:Tensión tangencial_2.png|900px|center]] 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
% Se realiza el mallado
h=0.1
r=1:h:2;
t=pi/4:h:3*pi/4;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
% Tensión tangencial
tension = ((sin(2.*tt-pi/2)./2.*rr).*(1-log(rr)));
% Tensión tangencial en 2D
subplot(1,2,1)
surf(xx,yy,tension)
axis equal
shading interp
view(2)
colorbar
title('Tensión tangencial en la dirección del eje e_ρ');
% Tensión tangencial en 2D
subplot(1,2,2)
surf(xx,yy,tension)
axis equal
shading interp
view(3)
colorbar
title('Tensión tangencial en 3D en la dirección del eje e_ρ');
}}

===Tensión de Von Mises===
La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. 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 teoría expone que un material dúctil comenzará a ceder en una ubicación cuando la tensión de von Mises sea igual al límite de tensión. En resumen 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. 
Esta viene determinada por 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> 
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).
Para este supuesto, la tensión de Von Mises representará la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
 
A continuación se muestra una representación gráfica de la tensión de Von Mises sobre la placa:
[[Archivo:Von mises22b.png|900px||centre|]]
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
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);
zz=zeros(size(xx));

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

%matriz de tensiones
sig = [];
vm = zeros(length(t), length(r));
%definimos las componentes de la matriz
for i = 1:length(t)
for j = 1:length(r)
sig(1,1) = M11(RR(i,j),TT(i,j));
sig(1,2) = M12(RR(i,j),TT(i,j));
sig(1,3) = 0;
sig(2,1) = M21(RR(i,j),TT(i,j));
sig(2,2) = M22(RR(i,j),TT(i,j));
sig(2,3) = 0;
sig(3,1) = 0;
sig(3,2) = 0;
sig(3,3) = M11(RR(i,j),TT(i,j));

[w,d] = eig(sig);%calculamos los autovalores
vm(i,j) = sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);%aplicamos la fórmula

end
end
subplot(1,3,2);
surf(xx,yy,vm)
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

m = vm(1:size(vm,1),1);
e = xx(1:length(xx),3);

surf(xx,yy,vm);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');

maxvm = max(m);
%hayamos el valor máximo
for k = 1:length(m)
if m(k) == maxvm
plot3(e(k),0,maxvm,'xr','markersize',10)
end
end

xlabel('Eje x')
zlabel('Eje z')
txt= ['valor max = ' num2str(maxvm)];
text(-1.85,maxvm+0.03,0.33,txt)
hold off

subplot(1,3,1)

surf(xx,yy,vm);
view(2)
axis([-3,3,-1,3])
colorbar;
title('Tensión de Von Mises');
}}
== Campo de Fuerzas ==
El campo de fuerzas, que actúa sobre la placa y que es el causante del desplazamiento observado, se calcula con la ecuación de Lamé:
 
A continuación se muestra una representación gráfica del campo de fuerzas:

[[Archivo:fuerza22b.png|centre|]]
El código empleado ha sido el siguiente: 
{{matlab|codigo=
clear
clear
u=1:0.1:2;
v=linspace(pi/4,3*pi/4,50);
[ro,theta]=meshgrid(u,v);
xx=ro.*cos(theta);
yy=ro.*sin(theta);
%la fuerza es -2laplaciano-gradiente(rotacional u)
fuerzai=-(((3.*cos(2.*theta).*log10(ro))./(2.*(ro.^2))).*sin(theta)-((log10((ro)+1).*sin(2.*theta))./(ro.^2)).*cos(theta))-2.*((-2.*log10(ro).*cos(2.*theta))./(ro.^2)).*sin(theta)
fuerzaj=-(((3.*cos(2.*theta).*log10(ro))./(2.*(ro.^2))).*cos(theta)-((log10((ro)+1).*sin(2.*theta))./(ro.^2)).*sin(theta))-2.*((2.*log10(ro).*cos(2.*theta))./(ro.^2)).*cos(theta)
quiver(xx,yy,fuerzai,fuerzaj);
}}

[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]

Visualización de campos escalares y vectoriales en elasticidad (Grupo 22-B)

2022-12-08T13:58:11Z

Sergio puente:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández}}
{{ Beta }}
== Introducción ==
El presente artículo tiene por objeto estudiar el desplazamiento experimentado por una pieza tras la aplicación de una fuerza.
 Para ello, se tienen dos cantidades físicas dependientes de las variables ρ y θ:
* La temperatura T(ρ,θ)
* El campo de desplazamientos <math>\vec u(ρ,θ)</math>, producido por la acción de una fuerza.
Para ello, se hará uso del software de programación y cálculo numérico Matlab/Octave.

== Descripción de la placa ==
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa de espesor diferencial con forma de cuarto de anillo circular, centrada en el origen y comprendida entre los radios 1 y 2 del plano <math>y≥|x|</math>, es decir, <math>1≤x^2+y^2≤4</math>.
Como se va a trabajar en coordenadas cilíndricas, esta región equivale a: <math>(ρ,θ)=[1,2]×[0,π/2]</math>. 
A continuación se muestra la sección de la placa.
[[Archivo:Placa plana.png|centre|]]
 
El código empleado para dibujar la placa es el siguiente:
{{matlab|codigo=

h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
mesh(x,y,z); %Dibujamos la placa
axis equal
axis([-3,3,-1,3]); %Ejes
view(2); %Vista en planta
title("Placa plana") %Nombre de la figura
}}

== Temperatura ==
La temperatura viene determinada por la siguiente función: <math> T(x,y)=log((x-3)^2)+2)</math>.
 En primer lugar se dibujarán las curvas de nivel de la temperatura.
[[Archivo:Curvas de nivel del campo de temperaturas.png|centre|]]
 El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel
axis equal;
axis([-3,3,-1,3]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
colorbar
}}
Ahora se dibujará la variación de la temperatura a lo largo de la placa.
[[Archivo:Campo de temperaturas.png|centre|]]
 '''El programa nos indica que la temperatura máxima alcanzada en la placa es de 3,0244ºC.'''
 El código empleado para dibujar el campo de temperatura a lo largo de la placa, así como para calcular la temperatura máxima alcanzada en la misma ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
surf(x,y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)
}}

== Gradiente de temperatura ==
Como se ha mencionado con anterioridad, el campo escalar que gobierna la temperatura en la placa es <math> T(x,y)=log((x-3)^2)+2)</math>.
 
El gradiente de un campo vectorial en coordenadas cartesianas se define como: 
<math>\nabla T(x,y,z) =\frac{d∂}{dx} +\frac{d∂}{dy} + \frac{d∂}{dz}</math>

De lo que resulta: <math>\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=>\frac{2(x-3)}{(x-3)^2+2} \vec i</math>
 
A continuación se muestra una representación del gradiente de la temperatura, como se puede apreciar este es ortogonal a las curvas de nivel de la temperatura ya que marca la dirección de máximo crecimiento de la misma.
Para este caso se observa que la temperatura crece en la dirección -<math>\vec i</math>
[[Archivo:Gradiente de temperatura_1.png|500px||center|]]
 
El código empleado para representar el gradiente ha sido el siguiente:
{{matlab|codigo=
h=0.2; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;;
T=log(((x-3).^2)+2); %Función de temperatura
grad=(2.*(x-3))./((x-3).^2+2); %Gradiente de temperatura
hold on
quiver(x,y,grad,0.*grad); %Se dibuja el gradiente
contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
axis([-3,3,-1,3]);
title("Gradiente de temperatura");
hold off
}}

== Campo de vectores ==
Se va a realizar una representación de el campo de desplazamientos que viene gobernado por la expresión: <math> \vec r(ρ,θ)= \frac{log(ρ)}{2}(sin(2θ-\frac{π}{2}))</math>
[[Archivo:Campo de desplazamientos sobre la curva.png|600px||center|]]
 
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=1/4*pi:pi/50:3/4*pi; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*(-sin(theta)); %Paso de la coordenada x del campo a cartesianas
FY=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*cos(theta); %Paso de la coordenada y del campo a cartesianas
quiver(X,Y,FX,FY); %Dibuja el campo vectorial
axis equal
axis([-3,3,-1,3]);
view(2) %Vista en planta
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')
}}
 

== Comparación de la placa antes y después del desplazamiento ==
En este apartado se va a realizar una comparación visual de la placa antes y después de sufrir el desplazamiento.
[[Archivo:Placa antes y después.png|800px||center|]]
[[Archivo:Comparación de placas.png|600px||center|]] 
El código empleado para llevar a cabo la comparación se expone a continuación:
{{matlab|codigo=
h=0.2;%paso de muestreo
r=1:h:2;
t=pi/4:h:((3*pi/4)+0.2);
[ro,theta]=meshgrid(r,t); %creamos el mallado
x=ro.*cos(theta);
y=ro.*sin(theta);
%sólido antes del desplazamieneto
subplot(1,2,1);
i=mesh(x,y,x.*0);
view(2);
set(i,'EdgeColor','r')
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa no desplazada');
%sólido despues del desplazamiento
subplot(1,2,2);
%pasamos el campo de desplazamiento a cartesianas
A=-(log10(ro)./2).*(sin(2.*theta-(pi/2)).*(sin(theta)));
B=(log10(ro)./2).*(cos(theta)).*(sin(2.*theta-(pi/2)));
%Este campo de desplazamiento no desplaza prácticamente la placa, si no que concentra la masa de la misma en el centroide.
X=x+A;
Y=y+B;
j=mesh(X,Y,0.*X);
view(2);
set(j,'EdgeColor','b');
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa desplazada');
%Comparación dde ambas placas
hold on
figure
j=mesh(X,Y,X*0);
axis equal
axis([-3,3,-1,3]);
view(2);
set(j,'EdgeColor','b');
title('Comparación de ambas placas: Roja (inicial), azul (final)');
hold on
i=mesh(x,y,x.*0);
set(i,'EdgeColor','r')
hold off
}}

==Divergencia==
La divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento. 
Para este supuesto se tiene que el campo de desplazamientos es <math> \vec u(ρ,θ) = \frac{ln(ρ)}{2}sin(2θ-\frac{π}{2})\vec e_θ </math> 
La divergencia resulta: <math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho·0})+\frac{\partial}{\partial{\theta}}(\frac{ln(ρ)}{2}sin(2θ-\frac{π}{2}))] </math>, es decir, <math>\nabla\cdot\vec u=\frac{ln(ρ)}{ρ}sin(2θ)</math>
 
La representación gráfica de la divergencia del campo de desplazamientos resulta: 
[[Archivo:divergencia_3.png|700px||center|]] 
El la placa, el valor máximo de la divergencia del campo de desplazamientos es 0.3466 y se alcanza en el punto (1.41421,1.41421). 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
%Paso de muestreo y definición de variables
h= 0.1;
r= 1:h:2;
t=pi/4:h:3*pi/4;
%Mallado
[rr,tt]= meshgrid(r,t);
%Parametrización
x=rr.*cos(tt);
y=rr.*sin(tt);
%Divergencia
div=(1./rr).*(log(rr).*sin(2.*(tt)));
%Divergencia en 2D
subplot(1,2,1)
surf(x,y,div)
shading interp
title('Divergencia en 2D')
axis([-3,3,-1,3]);
axis equal
colorbar;
view(2)
%Divergencia en 3D
subplot(1,2,2)
surf(x,y,div)
shading interp
axis([-3,3,-1,3]);
axis equal
colorbar;
view(3)
title('Divergencia en 3D')
maximo=max(max(div))
fprintf('El valor máximo de la divergencia es: %1.4f \n',maximo)
}}

== Rotacional ==
El rotacional de un campo vectorial permite observar la tendencia que tiene dicho campo a inducir rotación alrededor de un punto de la placa.
El rotacional, calculado en coordenadas cilíndricas viene determinado por la siguiente expresión:
 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|</math>. 
Si se aplica al campo de desplazamientos resulta: 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ·\frac{log(ρ}{2}·sin(2θ-\frac{π}{2}) & 0 \end{matrix}\right|=> -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}</math>. 
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
[[Archivo:Rotacional (2).png]] 
Como puede observarse, el rotacional es mayor en los puntos centrales de la placa, es por ello por lo que en la vista en 3 dimensiones la placa resulta doblada por esa zona. 
Del mismo modo, el programa nos indica que el valor máximo que alcanza el rotacional es 0.1504. 
El código empleado para dibujar el rotacional del campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
%Se crean los intervalos de las variables
h=50;
r=linspace(1,2,h);
t=linspace(pi/4,3*pi/4,h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=-(cos(2.*theta)./(2.*ro)).*log10((ro)+1)
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
shading interp
axis equal;
view(2);
colorbar;
title('Rotacional en 2D')
xlabel('Eje x');
ylabel('Eje y');
%Gráfica en 3D
subplot(1,2,2);
surf(x,y,rotacional);
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
Maximo=max(max(rotacional));
fprintf('el valor del rotacional máximo alcanzado es de: %1.4f \n',Maximo);
}}
==Tensor de deformaciones==
Las deformaciones que experimenta la placa se pueden modelizar mediante un tensor conocido como tensor de deformaciones <math>\epsilon(\vec u(ρ,θ))</math> que viene determinado por la parte simétrica del tensor <math> \nabla{\vec u(ρ,θ)} </math>, también llamado gradiente del campo de desplazamientos <math>\vec u(ρ,θ) </math>:
<math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> 
<math>\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}</math>
 
Lo primero será calcular <math>\nabla{\vec u(ρ,θ)}</math> que resulta: 
<math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial ρ} & \frac{\partial u_{1} }{\partial θ} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial ρ} & \frac{\partial u_{2} }{\partial θ} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial ρ} & \frac{\partial u_{3} }{\partial θ} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> 
 
 
::*<math>\frac{\partial \vec u}{\partial ρ}=(\frac{1}{2ρ})sen(2θ-\frac{π}{2}) \vec e_θ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial θ}= log(ρ)cos(2θ-\frac{π}{2}) \vec e_θ - (\frac{log(ρ)}{2})sen(2θ-\frac{π}{2}) \vec e_ρ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial z}= 0 \vec e_z</math>
 
 
El gradiente del campo de deformaciones resulta: 
<math>\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2})& 0\\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
A continuación se calcula la traspuesta de la matriz del gradiente del campo de deformaciones: 
<math>(\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
 
Aplicando la fórmula <math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> el tensor de deformaciones será igual a: 

<math>\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>

 
 

==Tensor de tensiones==
El tensor de tensiones viede deteminado por la siguiente expresión: <math>σ=λ·\nabla·\vec u I+2με</math> 
Desarrollando:
 
<math>\sigma = \left(\begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 & 0 \\ & \\0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) + 2 ·\left ( \begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix} \right) </math>
 
 
 
Finalmente resulta: <math>σ</math> = <math> \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right )
</math>
===Tensiones normales===
====Tensiones normales en la dirección del eje <math>\vec e_ρ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_ρ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math> 
 
Que resulta: <math>(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e rho.png|500px||left|]] [[Archivo:3drho_1.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
a=(log10(ro).*sin(2.*theta))./(ro);%elemento(1,1,1);
figure(1);
surf(x,y,a);
shading interp
axis equal
title('Tensión normal en la dirección e_ρ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_θ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_θ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math>
 
 
Que resulta: <math> (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e theta.png|500px|left]] [[Archivo:3drho.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
b=((log10(ro))./(ro))+2.*(log10(ro)).*sin(2.*theta);%elemento(2,2,2)
figure(1);
surf(x,y,b);
shading interp
axis equal
title('Tensión normal en la dirección e_θ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_z</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_z</math> es el siguiente: 
<math>|(σ·\vec e_θ))\vec e_θ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right ) ·\vec e_θ</math>
 
 
Que resulta: <math> (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
[[Archivo:Tension normal en e z.png|500px|left]] [[Archivo:3dz.png|500px|right]] 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha siso el siguiente: 
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
c=(log10(ro).*sin(2.*theta))./(ro);%elemento(3,3,3)
figure(1)
surf(x,y,c);
shading interp
axis equal
title('Tensión normal en la dirección e_k')
view(2);
colorbar;
}}

===Tensiones tangenciales respeto al plano ortogonal a <math>\vec e_ρ</math>===
Se van a calcular las tensiones tangenciales que sufre la placa en la dirección del eje <math>\vec e_ρ</math> en t=0. 
Para ello se conoce su el módulo es el siguiente: 
 
|<math>σ·\vec e_ρ-(\vec e_ρ·(σ·\vec e_ρ))\vec e_ρ ) </math>|. 
 
Desarrollando la expresión resulta: 
 
<math> \left | \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right )-
(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) \right | </math> 
 
Finalmente se obtiene: 
 
<math> \frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) </math>
La representación de la tensión tangencial que sufre la placa es la siguiente: 
[[Archivo:Tensión tangencial_2.png|900px|center]] 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
% Se realiza el mallado
h=0.1
r=1:h:2;
t=pi/4:h:3*pi/4;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
% Tensión tangencial
tension = ((sin(2.*tt-pi/2)./2.*rr).*(1-log(rr)));
% Tensión tangencial en 2D
subplot(1,2,1)
surf(xx,yy,tension)
axis equal
shading interp
view(2)
colorbar
title('Tensión tangencial en la dirección del eje e_ρ');
% Tensión tangencial en 2D
subplot(1,2,2)
surf(xx,yy,tension)
axis equal
shading interp
view(3)
colorbar
title('Tensión tangencial en 3D en la dirección del eje e_ρ');
}}

===Tensión de Von Mises===
La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. 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 teoría expone que un material dúctil comenzará a ceder en una ubicación cuando la tensión de von Mises sea igual al límite de tensión. En resumen 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. 
Esta viene determinada por 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> 
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).
Para este supuesto, la tensión de Von Mises representará la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
 
A continuación se muestra una representación gráfica de la tensión de Von Mises sobre la placa:
[[Archivo:Von mises22b.png|900px||centre|]]
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
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);
zz=zeros(size(xx));

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

%matriz de tensiones
sig = [];
vm = zeros(length(t), length(r));
%definimos las componentes de la matriz
for i = 1:length(t)
for j = 1:length(r)
sig(1,1) = M11(RR(i,j),TT(i,j));
sig(1,2) = M12(RR(i,j),TT(i,j));
sig(1,3) = 0;
sig(2,1) = M21(RR(i,j),TT(i,j));
sig(2,2) = M22(RR(i,j),TT(i,j));
sig(2,3) = 0;
sig(3,1) = 0;
sig(3,2) = 0;
sig(3,3) = M11(RR(i,j),TT(i,j));

[w,d] = eig(sig);%calculamos los autovalores
vm(i,j) = sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);%aplicamos la fórmula

end
end
subplot(1,3,2);
surf(xx,yy,vm)
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

m = vm(1:size(vm,1),1);
e = xx(1:length(xx),3);

surf(xx,yy,vm);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');

maxvm = max(m);
%hayamos el valor máximo
for k = 1:length(m)
if m(k) == maxvm
plot3(e(k),0,maxvm,'xr','markersize',10)
end
end

xlabel('Eje x')
zlabel('Eje z')
txt= ['valor max = ' num2str(maxvm)];
text(-1.85,maxvm+0.03,0.33,txt)
hold off

subplot(1,3,1)

surf(xx,yy,vm);
view(2)
axis([-3,3,-1,3])
colorbar;
title('Tensión de Von Mises');
}}
== Campo de Fuerzas ==
El campo de fuerzas, que actúa sobre la placa y que es el causante del desplazamiento observado, se calcula con la ecuación de Lamé:
 
A continuación se muestra una representación gráfica del campo de fuerzas:

[[Archivo:fuerza22b.png|centre|]]
El código empleado ha sido el siguiente: 
{{matlab|codigo=
clear
clear
u=1:0.1:2;
v=linspace(pi/4,3*pi/4,50);
[ro,theta]=meshgrid(u,v);
xx=ro.*cos(theta);
yy=ro.*sin(theta);
fuerzai=-(((3.*cos(2.*theta).*log10(ro))./(2.*(ro.^2))).*sin(theta)-((log10((ro)+1).*sin(2.*theta))./(ro.^2)).*cos(theta))-2.*((-2.*log10(ro).*cos(2.*theta))./(ro.^2)).*sin(theta)
fuerzaj=-(((3.*cos(2.*theta).*log10(ro))./(2.*(ro.^2))).*cos(theta)-((log10((ro)+1).*sin(2.*theta))./(ro.^2)).*sin(theta))-2.*((2.*log10(ro).*cos(2.*theta))./(ro.^2)).*cos(theta)
quiver(xx,yy,fuerzai,fuerzaj);
}}

[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]

Visualización de campos escalares y vectoriales en elasticidad (Grupo 22-B)

2022-12-08T13:56:27Z

Sergio puente:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández}}
{{ Beta }}
== Introducción ==
El presente artículo tiene por objeto estudiar el desplazamiento experimentado por una pieza tras la aplicación de una fuerza.
 Para ello, se tienen dos cantidades físicas dependientes de las variables ρ y θ:
* La temperatura T(ρ,θ)
* El campo de desplazamientos <math>\vec u(ρ,θ)</math>, producido por la acción de una fuerza.
Para ello, se hará uso del software de programación y cálculo numérico Matlab/Octave.

== Descripción de la placa ==
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa de espesor diferencial con forma de cuarto de anillo circular, centrada en el origen y comprendida entre los radios 1 y 2 del plano <math>y≥|x|</math>, es decir, <math>1≤x^2+y^2≤4</math>.
Como se va a trabajar en coordenadas cilíndricas, esta región equivale a: <math>(ρ,θ)=[1,2]×[0,π/2]</math>. 
A continuación se muestra la sección de la placa.
[[Archivo:Placa plana.png|centre|]]
 
El código empleado para dibujar la placa es el siguiente:
{{matlab|codigo=

h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
mesh(x,y,z); %Dibujamos la placa
axis equal
axis([-3,3,-1,3]); %Ejes
view(2); %Vista en planta
title("Placa plana") %Nombre de la figura
}}

== Temperatura ==
La temperatura viene determinada por la siguiente función: <math> T(x,y)=log((x-3)^2)+2)</math>.
 En primer lugar se dibujarán las curvas de nivel de la temperatura.
[[Archivo:Curvas de nivel del campo de temperaturas.png|centre|]]
 El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel
axis equal;
axis([-3,3,-1,3]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
colorbar
}}
Ahora se dibujará la variación de la temperatura a lo largo de la placa.
[[Archivo:Campo de temperaturas.png|centre|]]
 '''El programa nos indica que la temperatura máxima alcanzada en la placa es de 3,0244ºC.'''
 El código empleado para dibujar el campo de temperatura a lo largo de la placa, así como para calcular la temperatura máxima alcanzada en la misma ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
surf(x,y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)
}}

== Gradiente de temperatura ==
Como se ha mencionado con anterioridad, el campo escalar que gobierna la temperatura en la placa es <math> T(x,y)=log((x-3)^2)+2)</math>.
 
El gradiente de un campo vectorial en coordenadas cartesianas se define como: 
<math>\nabla T(x,y,z) =\frac{d∂}{dx} +\frac{d∂}{dy} + \frac{d∂}{dz}</math>

De lo que resulta: <math>\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=>\frac{2(x-3)}{(x-3)^2+2} \vec i</math>
 
A continuación se muestra una representación del gradiente de la temperatura, como se puede apreciar este es ortogonal a las curvas de nivel de la temperatura ya que marca la dirección de máximo crecimiento de la misma.
Para este caso se observa que la temperatura crece en la dirección -<math>\vec i</math>
[[Archivo:Gradiente de temperatura_1.png|500px||center|]]
 
El código empleado para representar el gradiente ha sido el siguiente:
{{matlab|codigo=
h=0.2; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;;
T=log(((x-3).^2)+2); %Función de temperatura
grad=(2.*(x-3))./((x-3).^2+2); %Gradiente de temperatura
hold on
quiver(x,y,grad,0.*grad); %Se dibuja el gradiente
contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
axis([-3,3,-1,3]);
title("Gradiente de temperatura");
hold off
}}

== Campo de vectores ==
Se va a realizar una representación de el campo de desplazamientos que viene gobernado por la expresión: <math> \vec r(ρ,θ)= \frac{log(ρ)}{2}(sin(2θ-\frac{π}{2}))</math>
[[Archivo:Campo de desplazamientos sobre la curva.png|600px||center|]]
 
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=1/4*pi:pi/50:3/4*pi; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*(-sin(theta)); %Paso de la coordenada x del campo a cartesianas
FY=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*cos(theta); %Paso de la coordenada y del campo a cartesianas
quiver(X,Y,FX,FY); %Dibuja el campo vectorial
axis equal
axis([-3,3,-1,3]);
view(2) %Vista en planta
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')
}}
 

== Comparación de la placa antes y después del desplazamiento ==
En este apartado se va a realizar una comparación visual de la placa antes y después de sufrir el desplazamiento.
[[Archivo:Placa antes y después.png|800px||center|]]
[[Archivo:Comparación de placas.png|600px||center|]] 
El código empleado para llevar a cabo la comparación se expone a continuación:
{{matlab|codigo=
h=0.2;%paso de muestreo
r=1:h:2;
t=pi/4:h:((3*pi/4)+0.2);
[ro,theta]=meshgrid(r,t); %creamos el mallado
x=ro.*cos(theta);
y=ro.*sin(theta);
%sólido antes del desplazamieneto
subplot(1,2,1);
i=mesh(x,y,x.*0);
view(2);
set(i,'EdgeColor','r')
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa no desplazada');
%sólido despues del desplazamiento
subplot(1,2,2);
%pasamos el campo de desplazamiento a cartesianas
A=-(log10(ro)./2).*(sin(2.*theta-(pi/2)).*(sin(theta)));
B=(log10(ro)./2).*(cos(theta)).*(sin(2.*theta-(pi/2)));
%Este campo de desplazamiento no desplaza prácticamente la placa, si no que concentra la masa de la misma en el centroide.
X=x+A;
Y=y+B;
j=mesh(X,Y,0.*X);
view(2);
set(j,'EdgeColor','b');
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa desplazada');
%Comparación dde ambas placas
hold on
figure
j=mesh(X,Y,X*0);
axis equal
axis([-3,3,-1,3]);
view(2);
set(j,'EdgeColor','b');
title('Comparación de ambas placas: Roja (inicial), azul (final)');
hold on
i=mesh(x,y,x.*0);
set(i,'EdgeColor','r')
hold off
}}

==Divergencia==
La divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento. 
Para este supuesto se tiene que el campo de desplazamientos es <math> \vec u(ρ,θ) = \frac{ln(ρ)}{2}sin(2θ-\frac{π}{2})\vec e_θ </math> 
La divergencia resulta: <math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho·0})+\frac{\partial}{\partial{\theta}}(\frac{ln(ρ)}{2}sin(2θ-\frac{π}{2}))] </math>, es decir, <math>\nabla\cdot\vec u=\frac{ln(ρ)}{ρ}sin(2θ)</math>
 
La representación gráfica de la divergencia del campo de desplazamientos resulta: 
[[Archivo:divergencia_3.png|700px||center|]] 
El la placa, el valor máximo de la divergencia del campo de desplazamientos es 0.3466 y se alcanza en el punto (1.41421,1.41421). 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
%Paso de muestreo y definición de variables
h= 0.1;
r= 1:h:2;
t=pi/4:h:3*pi/4;
%Mallado
[rr,tt]= meshgrid(r,t);
%Parametrización
x=rr.*cos(tt);
y=rr.*sin(tt);
%Divergencia
div=(1./rr).*(log(rr).*sin(2.*(tt)));
%Divergencia en 2D
subplot(1,2,1)
surf(x,y,div)
shading interp
title('Divergencia en 2D')
axis([-3,3,-1,3]);
axis equal
colorbar;
view(2)
%Divergencia en 3D
subplot(1,2,2)
surf(x,y,div)
shading interp
axis([-3,3,-1,3]);
axis equal
colorbar;
view(3)
title('Divergencia en 3D')
maximo=max(max(div))
fprintf('El valor máximo de la divergencia es: %1.4f \n',maximo)
}}

== Rotacional ==
El rotacional de un campo vectorial permite observar la tendencia que tiene dicho campo a inducir rotación alrededor de un punto de la placa.
El rotacional, calculado en coordenadas cilíndricas viene determinado por la siguiente expresión:
 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|</math>. 
Si se aplica al campo de desplazamientos resulta: 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ·\frac{log(ρ}{2}·sin(2θ-\frac{π}{2}) & 0 \end{matrix}\right|=> -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}</math>. 
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
[[Archivo:Rotacional (2).png]] 
Como puede observarse, el rotacional es mayor en los puntos centrales de la placa, es por ello por lo que en la vista en 3 dimensiones la placa resulta doblada por esa zona. 
Del mismo modo, el programa nos indica que el valor máximo que alcanza el rotacional es 0.1504. 
El código empleado para dibujar el rotacional del campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
%Se crean los intervalos de las variables
h=50;
r=linspace(1,2,h);
t=linspace(pi/4,3*pi/4,h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=-(cos(2.*theta)./(2.*ro)).*log10((ro)+1)
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
shading interp
axis equal;
view(2);
colorbar;
title('Rotacional en 2D')
xlabel('Eje x');
ylabel('Eje y');
%Gráfica en 3D
subplot(1,2,2);
surf(x,y,rotacional);
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
Maximo=max(max(rotacional));
fprintf('el valor del rotacional máximo alcanzado es de: %1.4f \n',Maximo);
}}
==Tensor de deformaciones==
Las deformaciones que experimenta la placa se pueden modelizar mediante un tensor conocido como tensor de deformaciones <math>\epsilon(\vec u(ρ,θ))</math> que viene determinado por la parte simétrica del tensor <math> \nabla{\vec u(ρ,θ)} </math>, también llamado gradiente del campo de desplazamientos <math>\vec u(ρ,θ) </math>:
<math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> 
<math>\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}</math>
 
Lo primero será calcular <math>\nabla{\vec u(ρ,θ)}</math> que resulta: 
<math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial ρ} & \frac{\partial u_{1} }{\partial θ} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial ρ} & \frac{\partial u_{2} }{\partial θ} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial ρ} & \frac{\partial u_{3} }{\partial θ} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> 
 
 
::*<math>\frac{\partial \vec u}{\partial ρ}=(\frac{1}{2ρ})sen(2θ-\frac{π}{2}) \vec e_θ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial θ}= log(ρ)cos(2θ-\frac{π}{2}) \vec e_θ - (\frac{log(ρ)}{2})sen(2θ-\frac{π}{2}) \vec e_ρ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial z}= 0 \vec e_z</math>
 
 
El gradiente del campo de deformaciones resulta: 
<math>\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2})& 0\\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
A continuación se calcula la traspuesta de la matriz del gradiente del campo de deformaciones: 
<math>(\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
 
Aplicando la fórmula <math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> el tensor de deformaciones será igual a: 

<math>\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>

 
 

==Tensor de tensiones==
El tensor de tensiones viede deteminado por la siguiente expresión: <math>σ=λ·\nabla·\vec u I+2με</math> 
Desarrollando:
 
<math>\sigma = \left(\begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 & 0 \\ & \\0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) + 2 ·\left ( \begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix} \right) </math>
 
 
 
Finalmente resulta: <math>σ</math> = <math> \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right )
</math>
===Tensiones normales===
====Tensiones normales en la dirección del eje <math>\vec e_ρ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_ρ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math> 
 
Que resulta: <math>(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e rho.png|500px||left|]] [[Archivo:3drho_1.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
a=(log10(ro).*sin(2.*theta))./(ro);%elemento(1,1,1);
figure(1);
surf(x,y,a);
shading interp
axis equal
title('Tensión normal en la dirección e_ρ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_θ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_θ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math>
 
 
Que resulta: <math> (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e theta.png|500px|left]] [[Archivo:3drho.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
b=((log10(ro))./(ro))+2.*(log10(ro)).*sin(2.*theta);%elemento(2,2,2)
figure(1);
surf(x,y,b);
shading interp
axis equal
title('Tensión normal en la dirección e_θ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_z</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_z</math> es el siguiente: 
<math>|(σ·\vec e_θ))\vec e_θ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right ) ·\vec e_θ</math>
 
 
Que resulta: <math> (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
[[Archivo:Tension normal en e z.png|500px|left]] [[Archivo:3dz.png|500px|right]] 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha siso el siguiente: 
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
c=(log10(ro).*sin(2.*theta))./(ro);%elemento(3,3,3)
figure(1)
surf(x,y,c);
shading interp
axis equal
title('Tensión normal en la dirección e_k')
view(2);
colorbar;
}}

===Tensiones tangenciales respeto al plano ortogonal a <math>\vec e_ρ</math>===
Se van a calcular las tensiones tangenciales que sufre la placa en la dirección del eje <math>\vec e_ρ</math> en t=0. 
Para ello se conoce su el módulo es el siguiente: 
 
|<math>σ·\vec e_ρ-(\vec e_ρ·(σ·\vec e_ρ))\vec e_ρ ) </math>|. 
 
Desarrollando la expresión resulta: 
 
<math> \left | \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right )-
(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) \right | </math> 
 
Finalmente se obtiene: 
 
<math> \frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) </math>
La representación de la tensión tangencial que sufre la placa es la siguiente: 
[[Archivo:Tensión tangencial_2.png|900px|center]] 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
% Se realiza el mallado
h=0.1
r=1:h:2;
t=pi/4:h:3*pi/4;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
% Tensión tangencial
tension = ((sin(2.*tt-pi/2)./2.*rr).*(1-log(rr)));
% Tensión tangencial en 2D
subplot(1,2,1)
surf(xx,yy,tension)
axis equal
shading interp
view(2)
colorbar
title('Tensión tangencial en la dirección del eje e_ρ');
% Tensión tangencial en 2D
subplot(1,2,2)
surf(xx,yy,tension)
axis equal
shading interp
view(3)
colorbar
title('Tensión tangencial en 3D en la dirección del eje e_ρ');
}}

===Tensión de Von Mises===
La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. 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 teoría expone que un material dúctil comenzará a ceder en una ubicación cuando la tensión de von Mises sea igual al límite de tensión. En resumen 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. 
Esta viene determinada por 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> 
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).
Para este supuesto, la tensión de Von Mises representará la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
 
A continuación se muestra una representación gráfica de la tensión de Von Mises sobre la placa:
[[Archivo:Von mises22b.png|900px||centre|]]
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
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);
zz=zeros(size(xx));

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

%matriz de tensiones
sig = [];
vm = zeros(length(t), length(r));
%definimos las componentes de la matriz
for i = 1:length(t)
for j = 1:length(r)
sig(1,1) = M11(RR(i,j),TT(i,j));
sig(1,2) = M12(RR(i,j),TT(i,j));
sig(1,3) = 0;
sig(2,1) = M21(RR(i,j),TT(i,j));
sig(2,2) = M22(RR(i,j),TT(i,j));
sig(2,3) = 0;
sig(3,1) = 0;
sig(3,2) = 0;
sig(3,3) = M11(RR(i,j),TT(i,j));

[w,d] = eig(sig);%calculamos los autovalores
vm(i,j) = sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);%aplicamos la fórmula

end
end
subplot(1,3,2);
surf(xx,yy,vm)
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

m = vm(1:size(vm,1),1);
e = xx(1:length(xx),3);

surf(xx,yy,vm);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');

maxvm = max(m);
%hayamos el valor máximo
for k = 1:length(m)
if m(k) == maxvm
plot3(e(k),0,maxvm,'xr','markersize',10)
end
end

xlabel('Eje x')
zlabel('Eje z')
txt= ['valor max = ' num2str(maxvm)];
text(-1.85,maxvm+0.03,0.33,txt)
hold off

subplot(1,3,1)

surf(xx,yy,vm);
view(2)
axis([-3,3,-1,3])
colorbar;
title('Tensión de Von Mises');
}}
== Campo de Fuerzas ==
El campo de fuerzas, que actúa sobre la placa y que es el causante del desplazamiento observado, se calcula con la ecuación de Lamé:
 
A continuación se muestra una representación gráfica del campo de fuerzas:

[[Archivo:fuerza22b.png|centre|]]

[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]

Visualización de campos escalares y vectoriales en elasticidad (Grupo 22-B)

2022-12-08T13:53:24Z

Sergio puente:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández}}
{{ Beta }}
== Introducción ==
El presente artículo tiene por objeto estudiar el desplazamiento experimentado por una pieza tras la aplicación de una fuerza.
 Para ello, se tienen dos cantidades físicas dependientes de las variables ρ y θ:
* La temperatura T(ρ,θ)
* El campo de desplazamientos <math>\vec u(ρ,θ)</math>, producido por la acción de una fuerza.
Para ello, se hará uso del software de programación y cálculo numérico Matlab/Octave.

== Descripción de la placa ==
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa de espesor diferencial con forma de cuarto de anillo circular, centrada en el origen y comprendida entre los radios 1 y 2 del plano <math>y≥|x|</math>, es decir, <math>1≤x^2+y^2≤4</math>.
Como se va a trabajar en coordenadas cilíndricas, esta región equivale a: <math>(ρ,θ)=[1,2]×[0,π/2]</math>. 
A continuación se muestra la sección de la placa.
[[Archivo:Placa plana.png|centre|]]
 
El código empleado para dibujar la placa es el siguiente:
{{matlab|codigo=

h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
mesh(x,y,z); %Dibujamos la placa
axis equal
axis([-3,3,-1,3]); %Ejes
view(2); %Vista en planta
title("Placa plana") %Nombre de la figura
}}

== Temperatura ==
La temperatura viene determinada por la siguiente función: <math> T(x,y)=log((x-3)^2)+2)</math>.
 En primer lugar se dibujarán las curvas de nivel de la temperatura.
[[Archivo:Curvas de nivel del campo de temperaturas.png|centre|]]
 El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel
axis equal;
axis([-3,3,-1,3]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
colorbar
}}
Ahora se dibujará la variación de la temperatura a lo largo de la placa.
[[Archivo:Campo de temperaturas.png|centre|]]
 '''El programa nos indica que la temperatura máxima alcanzada en la placa es de 3,0244ºC.'''
 El código empleado para dibujar el campo de temperatura a lo largo de la placa, así como para calcular la temperatura máxima alcanzada en la misma ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
surf(x,y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)
}}

== Gradiente de temperatura ==
Como se ha mencionado con anterioridad, el campo escalar que gobierna la temperatura en la placa es <math> T(x,y)=log((x-3)^2)+2)</math>.
 
El gradiente de un campo vectorial en coordenadas cartesianas se define como: 
<math>\nabla T(x,y,z) =\frac{d∂}{dx} +\frac{d∂}{dy} + \frac{d∂}{dz}</math>

De lo que resulta: <math>\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=>\frac{2(x-3)}{(x-3)^2+2} \vec i</math>
 
A continuación se muestra una representación del gradiente de la temperatura, como se puede apreciar este es ortogonal a las curvas de nivel de la temperatura ya que marca la dirección de máximo crecimiento de la misma.
Para este caso se observa que la temperatura crece en la dirección -<math>\vec i</math>
[[Archivo:Gradiente de temperatura_1.png|500px||center|]]
 
El código empleado para representar el gradiente ha sido el siguiente:
{{matlab|codigo=
h=0.2; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;;
T=log(((x-3).^2)+2); %Función de temperatura
grad=(2.*(x-3))./((x-3).^2+2); %Gradiente de temperatura
hold on
quiver(x,y,grad,0.*grad); %Se dibuja el gradiente
contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
axis([-3,3,-1,3]);
title("Gradiente de temperatura");
hold off
}}

== Campo de vectores ==
Se va a realizar una representación de el campo de desplazamientos que viene gobernado por la expresión: <math> \vec r(ρ,θ)= \frac{log(ρ)}{2}(sin(2θ-\frac{π}{2}))</math>
[[Archivo:Campo de desplazamientos sobre la curva.png|600px||center|]]
 
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=1/4*pi:pi/50:3/4*pi; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*(-sin(theta)); %Paso de la coordenada x del campo a cartesianas
FY=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*cos(theta); %Paso de la coordenada y del campo a cartesianas
quiver(X,Y,FX,FY); %Dibuja el campo vectorial
axis equal
axis([-3,3,-1,3]);
view(2) %Vista en planta
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')
}}
 

== Comparación de la placa antes y después del desplazamiento ==
En este apartado se va a realizar una comparación visual de la placa antes y después de sufrir el desplazamiento.
[[Archivo:Placa antes y después.png|800px||center|]]
[[Archivo:Comparación de placas.png|600px||center|]] 
El código empleado para llevar a cabo la comparación se expone a continuación:
{{matlab|codigo=
h=0.2;%paso de muestreo
r=1:h:2;
t=pi/4:h:((3*pi/4)+0.2);
[ro,theta]=meshgrid(r,t); %creamos el mallado
x=ro.*cos(theta);
y=ro.*sin(theta);
%sólido antes del desplazamieneto
subplot(1,2,1);
i=mesh(x,y,x.*0);
view(2);
set(i,'EdgeColor','r')
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa no desplazada');
%sólido despues del desplazamiento
subplot(1,2,2);
%pasamos el campo de desplazamiento a cartesianas
A=-(log10(ro)./2).*(sin(2.*theta-(pi/2)).*(sin(theta)));
B=(log10(ro)./2).*(cos(theta)).*(sin(2.*theta-(pi/2)));
%Este campo de desplazamiento no desplaza prácticamente la placa, si no que concentra la masa de la misma en el centroide.
X=x+A;
Y=y+B;
j=mesh(X,Y,0.*X);
view(2);
set(j,'EdgeColor','b');
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa desplazada');
%Comparación dde ambas placas
hold on
figure
j=mesh(X,Y,X*0);
axis equal
axis([-3,3,-1,3]);
view(2);
set(j,'EdgeColor','b');
title('Comparación de ambas placas: Roja (inicial), azul (final)');
hold on
i=mesh(x,y,x.*0);
set(i,'EdgeColor','r')
hold off
}}

==Divergencia==
La divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento. 
Para este supuesto se tiene que el campo de desplazamientos es <math> \vec u(ρ,θ) = \frac{ln(ρ)}{2}sin(2θ-\frac{π}{2})\vec e_θ </math> 
La divergencia resulta: <math>\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho·0})+\frac{\partial}{\partial{\theta}}(\frac{ln(ρ)}{2}sin(2θ-\frac{π}{2}))] </math>, es decir, <math>\nabla\cdot\vec u=\frac{ln(ρ)}{ρ}sin(2θ)</math>
 
La representación gráfica de la divergencia del campo de desplazamientos resulta: 
[[Archivo:divergencia_3.png|700px||center|]] 
El la placa, el valor máximo de la divergencia del campo de desplazamientos es 0.3466 y se alcanza en el punto (1.41421,1.41421). 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
%Paso de muestreo y definición de variables
h= 0.1;
r= 1:h:2;
t=pi/4:h:3*pi/4;
%Mallado
[rr,tt]= meshgrid(r,t);
%Parametrización
x=rr.*cos(tt);
y=rr.*sin(tt);
%Divergencia
div=(1./rr).*(log(rr).*sin(2.*(tt)));
%Divergencia en 2D
subplot(1,2,1)
surf(x,y,div)
shading interp
title('Divergencia en 2D')
axis([-3,3,-1,3]);
axis equal
colorbar;
view(2)
%Divergencia en 3D
subplot(1,2,2)
surf(x,y,div)
shading interp
axis([-3,3,-1,3]);
axis equal
colorbar;
view(3)
title('Divergencia en 3D')
maximo=max(max(div))
fprintf('El valor máximo de la divergencia es: %1.4f \n',maximo)
}}

== Rotacional ==
El rotacional de un campo vectorial permite observar la tendencia que tiene dicho campo a inducir rotación alrededor de un punto de la placa.
El rotacional, calculado en coordenadas cilíndricas viene determinado por la siguiente expresión:
 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|</math>. 
Si se aplica al campo de desplazamientos resulta: 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ·\frac{log(ρ}{2}·sin(2θ-\frac{π}{2}) & 0 \end{matrix}\right|=> -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}</math>. 
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
[[Archivo:Rotacional (2).png]] 
Como puede observarse, el rotacional es mayor en los puntos centrales de la placa, es por ello por lo que en la vista en 3 dimensiones la placa resulta doblada por esa zona. 
Del mismo modo, el programa nos indica que el valor máximo que alcanza el rotacional es 0.1504. 
El código empleado para dibujar el rotacional del campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
%Se crean los intervalos de las variables
h=50;
r=linspace(1,2,h);
t=linspace(pi/4,3*pi/4,h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=-(cos(2.*theta)./(2.*ro)).*log10((ro)+1)
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
shading interp
axis equal;
view(2);
colorbar;
title('Rotacional en 2D')
xlabel('Eje x');
ylabel('Eje y');
%Gráfica en 3D
subplot(1,2,2);
surf(x,y,rotacional);
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
Maximo=max(max(rotacional));
fprintf('el valor del rotacional máximo alcanzado es de: %1.4f \n',Maximo);
}}
==Tensor de deformaciones==
Las deformaciones que experimenta la placa se pueden modelizar mediante un tensor conocido como tensor de deformaciones <math>\epsilon(\vec u(ρ,θ))</math> que viene determinado por la parte simétrica del tensor <math> \nabla{\vec u(ρ,θ)} </math>, también llamado gradiente del campo de desplazamientos <math>\vec u(ρ,θ) </math>:
<math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> 
<math>\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}</math>
 
Lo primero será calcular <math>\nabla{\vec u(ρ,θ)}</math> que resulta: 
<math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial ρ} & \frac{\partial u_{1} }{\partial θ} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial ρ} & \frac{\partial u_{2} }{\partial θ} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial ρ} & \frac{\partial u_{3} }{\partial θ} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> 
 
 
::*<math>\frac{\partial \vec u}{\partial ρ}=(\frac{1}{2ρ})sen(2θ-\frac{π}{2}) \vec e_θ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial θ}= log(ρ)cos(2θ-\frac{π}{2}) \vec e_θ - (\frac{log(ρ)}{2})sen(2θ-\frac{π}{2}) \vec e_ρ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial z}= 0 \vec e_z</math>
 
 
El gradiente del campo de deformaciones resulta: 
<math>\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2})& 0\\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
A continuación se calcula la traspuesta de la matriz del gradiente del campo de deformaciones: 
<math>(\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
 
Aplicando la fórmula <math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> el tensor de deformaciones será igual a: 

<math>\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>

 
 

==Tensor de tensiones==
El tensor de tensiones viede deteminado por la siguiente expresión: <math>σ=λ·\nabla·\vec u I+2με</math> 
Desarrollando:
 
<math>\sigma = \left(\begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 & 0 \\ & \\0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) + 2 ·\left ( \begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix} \right) </math>
 
 
 
Finalmente resulta: <math>σ</math> = <math> \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right )
</math>
===Tensiones normales===
====Tensiones normales en la dirección del eje <math>\vec e_ρ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_ρ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math> 
 
Que resulta: <math>(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e rho.png|500px||left|]] [[Archivo:3drho_1.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
a=(log10(ro).*sin(2.*theta))./(ro);%elemento(1,1,1);
figure(1);
surf(x,y,a);
shading interp
axis equal
title('Tensión normal en la dirección e_ρ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_θ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_θ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math>
 
 
Que resulta: <math> (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e theta.png|500px|left]] [[Archivo:3drho.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
b=((log10(ro))./(ro))+2.*(log10(ro)).*sin(2.*theta);%elemento(2,2,2)
figure(1);
surf(x,y,b);
shading interp
axis equal
title('Tensión normal en la dirección e_θ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_z</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_z</math> es el siguiente: 
<math>|(σ·\vec e_θ))\vec e_θ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right ) ·\vec e_θ</math>
 
 
Que resulta: <math> (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
[[Archivo:Tension normal en e z.png|500px|left]] [[Archivo:3dz.png|500px|right]] 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha siso el siguiente: 
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
c=(log10(ro).*sin(2.*theta))./(ro);%elemento(3,3,3)
figure(1)
surf(x,y,c);
shading interp
axis equal
title('Tensión normal en la dirección e_k')
view(2);
colorbar;
}}

===Tensiones tangenciales respeto al plano ortogonal a <math>\vec e_ρ</math>===
Se van a calcular las tensiones tangenciales que sufre la placa en la dirección del eje <math>\vec e_ρ</math> en t=0. 
Para ello se conoce su el módulo es el siguiente: 
 
|<math>σ·\vec e_ρ-(\vec e_ρ·(σ·\vec e_ρ))\vec e_ρ ) </math>|. 
 
Desarrollando la expresión resulta: 
 
<math> \left | \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right )-
(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) \right | </math> 
 
Finalmente se obtiene: 
 
<math> \frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) </math>
La representación de la tensión tangencial que sufre la placa es la siguiente: 
[[Archivo:Tensión tangencial_2.png|900px|center]] 
El código empleado ha sido el siguiente: 
{{matlab|codigo=
% Se realiza el mallado
h=0.1
r=1:h:2;
t=pi/4:h:3*pi/4;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
% Tensión tangencial
tension = ((sin(2.*tt-pi/2)./2.*rr).*(1-log(rr)));
% Tensión tangencial en 2D
subplot(1,2,1)
surf(xx,yy,tension)
axis equal
shading interp
view(2)
colorbar
title('Tensión tangencial en la dirección del eje e_ρ');
% Tensión tangencial en 2D
subplot(1,2,2)
surf(xx,yy,tension)
axis equal
shading interp
view(3)
colorbar
title('Tensión tangencial en 3D en la dirección del eje e_ρ');
}}

===Tensión de Von Mises===
La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. 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 teoría expone que un material dúctil comenzará a ceder en una ubicación cuando la tensión de von Mises sea igual al límite de tensión. En resumen 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. 
Esta viene determinada por 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> 
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).
Para este supuesto, la tensión de Von Mises representará la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
 
A continuación se muestra una representación gráfica de la tensión de Von Mises sobre la placa:
[[Archivo:Von mises22b.png|900px||centre|]]
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
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);
zz=zeros(size(xx));

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

%matriz de tensiones
sig = [];
vm = zeros(length(t), length(r));
%definimos las componentes de la matriz
for i = 1:length(t)
for j = 1:length(r)
sig(1,1) = M11(RR(i,j),TT(i,j));
sig(1,2) = M12(RR(i,j),TT(i,j));
sig(1,3) = 0;
sig(2,1) = M21(RR(i,j),TT(i,j));
sig(2,2) = M22(RR(i,j),TT(i,j));
sig(2,3) = 0;
sig(3,1) = 0;
sig(3,2) = 0;
sig(3,3) = M11(RR(i,j),TT(i,j));

[w,d] = eig(sig);%calculamos los autovalores
vm(i,j) = sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);%aplicamos la fórmula

end
end
subplot(1,3,2);
surf(xx,yy,vm)
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

m = vm(1:size(vm,1),1);
e = xx(1:length(xx),3);

surf(xx,yy,vm);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');

maxvm = max(m);
%hayamos el valor máximo
for k = 1:length(m)
if m(k) == maxvm
plot3(e(k),0,maxvm,'xr','markersize',10)
end
end

xlabel('Eje x')
zlabel('Eje z')
txt= ['valor max = ' num2str(maxvm)];
text(-1.85,maxvm+0.03,0.33,txt)
hold off

subplot(1,3,1)

surf(xx,yy,vm);
view(2)
axis([-3,3,-1,3])
colorbar;
title('Tensión de Von Mises');
}}
== Campo de Fuerzas ==
El campo de fuerzas, que actúa sobre la placa y que es el causante del desplazamiento observado, se calcula con la ecuación de Lamé:
 
A continuación se muestra una representación gráfica del campo de fuerzas:
[[Archivo:fuerza22b||centre|]]

[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Teoría de Campos]]
[[Categoría:TC22/23]]

Archivo:Fuerza22b.png

2022-12-08T13:52:42Z

Sergio puente:

Visualización de campos escalares y vectoriales en elasticidad (Grupo 22-B)

2022-12-08T08:17:03Z

Sergio puente:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández}}
{{ Beta }}
== Introducción ==
El presente artículo tiene por objeto estudiar el desplazamiento experimentado por una pieza tras la aplicación de una fuerza.
 Para ello, se tienen dos cantidades físicas dependientes de las variables ρ y θ:
* La temperatura T(ρ,θ)
* El campo de desplazamientos <math>\vec u(ρ,θ)</math>, producido por la acción de una fuerza.
Para ello, se hará uso del software de programación y cálculo numérico Matlab/Octave.

== Descripción de la placa ==
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa de espesor diferencial con forma de cuarto de anillo circular, centrada en el origen y comprendida entre los radios 1 y 2 del plano <math>y≥|x|</math>, es decir, <math>1≤x^2+y^2≤4</math>.
Como se va a trabajar en coordenadas cilíndricas, esta región equivale a: <math>(ρ,θ)=[1,2]×[0,π/2]</math>.
A continuación se muestra la sección de la placa.
[[Archivo:Placa plana.png|centre|]]
 
El código empleado para dibujar la placa es el siguiente:
{{matlab|codigo=

h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
mesh(x,y,z); %Dibujamos la placa
axis equal
axis([-3,3,-1,3]); %Ejes
view(2); %Vista en planta
title("Placa plana") %Nombre de la figura
}}

== Temperatura ==
La temperatura viene determinada por la siguiente función: <math> T(x,y)=log((x-3)^2)+2)</math>.
 En primer lugar se dibujarán las curvas de nivel de la temperatura.
[[Archivo:Curvas de nivel del campo de temperaturas.png|centre|]]
 El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel
axis equal;
axis([-3,3,-1,3]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
colorbar
}}
Ahora se dibujará la variación de la temperatura a lo largo de la placa.
[[Archivo:Campo de temperaturas.png|centre|]]
 '''El programa nos indica que la temperatura máxima alcanzada en la placa es de 3,0244ºC.'''
 El código empleado para dibujar el campo de temperatura a lo largo de la placa, así como para calcular la temperatura máxima alcanzada en la misma ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
surf(x,y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)
}}

== Gradiente de temperatura ==
Como se ha mencionado con anterioridad, el campo escalar que gobierna la temperatura en la placa es <math> T(x,y)=log((x-3)^2)+2)</math>.
 
El gradiente de un campo vectorial se define como: 
<center><math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial x} & \frac{\partial u_{1} }{\partial y} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial x} & \frac{\partial u_{2} }{\partial y} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial x} & \frac{\partial u_{3} }{\partial y} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> </center> 
Al no tener el campo componentes en y ni en z, la matriz resulta: 
<center> <math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial x} & 0 & 0\\
\frac{\partial u_{2} }{\partial x} & 0 & 0 \\
\frac{\partial u_{3} }{\partial x} & 0 & 0
\end{pmatrix} </math> </center>

De lo que resulta: <math>\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=>\frac{2(x-3)}{(x-3)^2+2} \vec i</math>
 
A continuación se muestra una representación del gradiente de la temperatura, como se puede apreciar este es ortogonal a las curvas de nivel de la temperatura ya que marca la dirección de máximo crecimiento de la misma.
Para este caso se observa que la temperatura crece en la dirección -<math>\vec i</math>
[[Archivo:Gradiente de temperatura_1.png]]
 
El código empleado para representar el gradiente ha sido el siguiente:
{{matlab|codigo=
h=0.2; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;;
T=log(((x-3).^2)+2); %Función de temperatura
grad=(2.*(x-3))./((x-3).^2+2); %Gradiente de temperatura
hold on
quiver(x,y,grad,0.*grad); %Se dibuja el gradiente
contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
axis([-3,3,-1,3]);
title("Gradiente de temperatura");
hold off
}}

== Campo de vectores ==
Se va a realizar una representación de el campo de desplazamientos que viene gobernado por la expresión: <math> \vec r(ρ,θ)= \frac{log(ρ)}{2}(sin(2θ-\frac{π}{2}))</math>
[[Archivo:Campo de desplazamientos sobre la curva.png]]
 
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=1/4*pi:pi/50:3/4*pi; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*(-sin(theta)); %Paso de la coordenada x del campo a cartesianas
FY=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*cos(theta); %Paso de la coordenada y del campo a cartesianas
quiver(X,Y,FX,FY); %Dibuja el campo vectorial
axis equal
axis([-3,3,-1,3]);
view(2) %Vista en planta
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')
}}
 
== Comparación de la placa antes y después del desplazamiento ==
En este apartado se va a realizar una comparación visual de la placa antes y después de sufrir el desplazamiento.
[[Archivo:Placa antes y después.png]]
[[Archivo:Comparación de placas.png]] 
{{matlab|codigo=
El código empleado para llevar a cabo la comparación se expone a continuación:
h=0.2;%paso de muestreo
r=1:h:2;
t=pi/4:h:((3*pi/4)+0.2);
[ro,theta]=meshgrid(r,t); %creamos el mallado
x=ro.*cos(theta);
y=ro.*sin(theta);
%sólido antes del desplazamieneto
subplot(1,2,1);
i=mesh(x,y,x.*0);
view(2);
set(i,'EdgeColor','r')
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa no desplazada');
%sólido despues del desplazamiento
subplot(1,2,2);
%pasamos el campo de desplazamiento a cartesianas
A=-(log10(ro)./2).*(sin(2.*theta-(pi/2)).*(sin(theta)));
B=(log10(ro)./2).*(cos(theta)).*(sin(2.*theta-(pi/2)));
%Este campo de desplazamiento no desplaza prácticamente la placa, si no que concentra la masa de la misma en el centroide.
X=x+A;
Y=y+B;
j=mesh(X,Y,0.*X);
view(2);
set(j,'EdgeColor','b');
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa desplazada');
%Comparación dde ambas placas
hold on
figure
j=mesh(X,Y,X*0);
axis equal
axis([-3,3,-1,3]);
view(2);
set(j,'EdgeColor','b');
title('Comparación de ambas placas: Roja (inicial), azul (final)');
hold on
i=mesh(x,y,x.*0);
set(i,'EdgeColor','r')
hold off
}}
== Rotacional ==
El rotacional de un campo vectorial permite observar la tendencia que tiene dicho campo a inducir rotación alrededor de un punto de la placa.
El rotacional, calculado en coordenadas cilíndricas viene determinado por la siguiente expresión:
 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|</math>. 
Si se aplica al campo de desplazamientos resulta: 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ·\frac{log(ρ}{2}·sin(2θ-\frac{π}{2}) & 0 \end{matrix}\right|=> -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}</math>. 
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
[[Archivo:Rotacional (2).png]] 
Como puede observarse, el rotacional es mayor en los puntos centrales de la placa, es por ello por lo que en la vista en 3 dimensiones la placa resulta doblada por esa zona. 
Del mismo modo, el programa nos indica que el valor máximo que alcanza el rotacional es 0.1504. 
El código empleado para dibujar el rotacional del campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
%Se crean los intervalos de las variables
h=50;
r=linspace(1,2,h);
t=linspace(pi/4,3*pi/4,h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=-(cos(2.*theta)./(2.*ro)).*log10((ro)+1)
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
shading interp
axis equal;
view(2);
colorbar;
title('Rotacional en 2D')
xlabel('Eje x');
ylabel('Eje y');
%Gráfica en 3D
subplot(1,2,2);
surf(x,y,rotacional);
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
Maximo=max(max(rotacional));
fprintf('el valor del rotacional máximo alcanzado es de: %1.4f \n',Maximo);
}}

==Tensor de deformaciones==
A continuación, se va a analizar el comportamiento de la placa al estar sometida a una serie de tensiones en un medio elástico lineal, isótropo y homogéneo. Se va a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento ondulatorio de los puntos de la misma dado por el siguiente campo vectorial: 
<math>\vec u(x,y,t)=\vec asin((k(\vec d \vec r -vt))</math>
 
Donde los desplazamientos permiten escribir el tensor de tensiones como: 
<math>σ=λ·\nabla·\vec u I+2με</math> 
Para ello, se va a definir 
<math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> 
Como la parte simétrica del tensor gradiente <math>\vec u</math> conocido como tensor de deformaciones. 
λ y μ son los llamados coeficientes de Lamé que por las características del material se conoce que son igual a 1.
Se define el tensor de deformaciones <math>\epsilon(\vec u(ρ,θ))</math> como la parte simétrica del tensor <math> \nabla{\vec u(ρ,θ)} </math>, llamado también gradiente del campo de desplazamientos <math>\vec u(ρ,θ) </math>:
: <math>\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}</math>.
 
Lo primero será calcular <math>\nabla{\vec u(ρ,θ)}</math> que resulta: 
<math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial ρ} & \frac{\partial u_{1} }{\partial θ} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial ρ} & \frac{\partial u_{2} }{\partial θ} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial ρ} & \frac{\partial u_{3} }{\partial θ} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> 
 
 
::*<math>\frac{\partial \vec u}{\partial ρ}=(\frac{1}{2ρ})sen(2θ-\frac{π}{2}) \vec e_θ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial θ}= log(ρ)cos(2θ-\frac{π}{2}) \vec e_θ - (\frac{log(ρ)}{2})sen(2θ-\frac{π}{2}) \vec e_ρ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial z}= 0 \vec e_z</math>
 
 
El gradiente del campo de deformaciones resulta: 
<math>\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2})& 0\\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
A continuación se calcula la traspuesta de la matriz del gradiente del campo de deformaciones: 
<math>(\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
 
Aplicando la fórmula <math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> el tensor de deformaciones será igual a: 

<math>\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>

 
 

===Tensiones normales===
====Tensiones normales en la dirección del eje <math>\vec e_ρ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_ρ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math> 
 
Que resulta: <math>(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e rho.png|500px||left|]] [[Archivo:3drho_1.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
a=(log10(ro).*sin(2.*theta))./(ro);%elemento(1,1,1);
figure(1);
surf(x,y,a);
shading interp
axis equal
title('Tensión normal en la dirección e_ρ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_θ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_θ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math>
 
 
Que resulta: <math> (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e theta.png|500px|left]] [[Archivo:3drho.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
b=((log10(ro))./(ro))+2.*(log10(ro)).*sin(2.*theta);%elemento(2,2,2)
figure(1);
surf(x,y,b);
shading interp
axis equal
title('Tensión normal en la dirección e_θ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_z</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_z</math> es el siguiente:
<math>|(σ·\vec e_θ))\vec e_θ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right ) ·\vec e_θ</math>
 
 
Que resulta: <math> (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
[[Archivo:Tension normal en e z.png|500px|left]] [[Archivo:3dz.png|500px|right]] 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha siso el siguiente: 
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
c=(log10(ro).*sin(2.*theta))./(ro);%elemento(3,3,3)
figure(1)
surf(x,y,c);
shading interp
axis equal
title('Tensión normal en la dirección e_k')
view(2);
colorbar;
}}

===Tensiones tangenciales respeto al plano ortogonal a <math>\vec e_ρ</math>===
===Tensión de Von Mises===
La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. 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 teoría expone que un material dúctil comenzará a ceder en una ubicación cuando la tensión de von Mises sea igual al límite de tensión. En resumen 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. 
Esta viene determinada por 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> 
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).
Para este supuesto, la tensión de Von Mises representará la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
 
A continuación se muestra la representación de dicha tensión:
[[Archivo:Von mises22b.png|centre|]]
 
El código empleado para la misma es:
{{matlab|codigo=
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);
zz=zeros(size(xx));

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

%matriz de tensiones
sig = [];
vm = zeros(length(t), length(r));
%definimos las componentes de la matriz
for i = 1:length(t)
for j = 1:length(r)
sig(1,1) = M11(RR(i,j),TT(i,j));
sig(1,2) = M12(RR(i,j),TT(i,j));
sig(1,3) = 0;
sig(2,1) = M21(RR(i,j),TT(i,j));
sig(2,2) = M22(RR(i,j),TT(i,j));
sig(2,3) = 0;
sig(3,1) = 0;
sig(3,2) = 0;
sig(3,3) = M11(RR(i,j),TT(i,j));

[w,d] = eig(sig);%calculamos los autovalores
vm(i,j) = sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);%aplicamos la fórmula

end
end
subplot(1,3,2);
surf(xx,yy,vm)
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

m = vm(1:size(vm,1),1);
e = xx(1:length(xx),3);

surf(xx,yy,vm);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');

maxvm = max(m);
%hayamos el valor máximo
for k = 1:length(m)
if m(k) == maxvm
plot3(e(k),0,maxvm,'xr','markersize',10)
end
end

xlabel('Eje x')
zlabel('Eje z')
txt= ['valor max = ' num2str(maxvm)];
text(-1.85,maxvm+0.03,0.33,txt)
hold off

subplot(1,3,1)

surf(xx,yy,vm);
view(2)
axis([-3,3,-1,3])
colorbar;
title('Tensión de Von Mises');
}}

[[Categoría:Teoría de Campos]]

Visualización de campos escalares y vectoriales en elasticidad (Grupo 22-B)

2022-12-08T08:15:27Z

Sergio puente:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC22/23|2022-23]] | Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández}}
{{ Beta }}
== Introducción ==
El presente artículo tiene por objeto estudiar el desplazamiento experimentado por una pieza tras la aplicación de una fuerza.
 Para ello, se tienen dos cantidades físicas dependientes de las variables ρ y θ:
* La temperatura T(ρ,θ)
* El campo de desplazamientos <math>\vec u(ρ,θ)</math>, producido por la acción de una fuerza.
Para ello, se hará uso del software de programación y cálculo numérico Matlab/Octave.

== Descripción de la placa ==
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa de espesor diferencial con forma de cuarto de anillo circular, centrada en el origen y comprendida entre los radios 1 y 2 del plano <math>y≥|x|</math>, es decir, <math>1≤x^2+y^2≤4</math>.
Como se va a trabajar en coordenadas cilíndricas, esta región equivale a: <math>(ρ,θ)=[1,2]×[0,π/2]</math>.
A continuación se muestra la sección de la placa.
[[Archivo:Placa plana.png|centre|]]
 
El código empleado para dibujar la placa es el siguiente:
{{matlab|codigo=

h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
mesh(x,y,z); %Dibujamos la placa
axis equal
axis([-3,3,-1,3]); %Ejes
view(2); %Vista en planta
title("Placa plana") %Nombre de la figura
}}

== Temperatura ==
La temperatura viene determinada por la siguiente función: <math> T(x,y)=log((x-3)^2)+2)</math>.
 En primer lugar se dibujarán las curvas de nivel de la temperatura.
[[Archivo:Curvas de nivel del campo de temperaturas.png|centre|]]
 El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel
axis equal;
axis([-3,3,-1,3]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
colorbar
}}
Ahora se dibujará la variación de la temperatura a lo largo de la placa.
[[Archivo:Campo de temperaturas.png|centre|]]
 '''El programa nos indica que la temperatura máxima alcanzada en la placa es de 3,0244ºC.'''
 El código empleado para dibujar el campo de temperatura a lo largo de la placa, así como para calcular la temperatura máxima alcanzada en la misma ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; % Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;
T=log(((x-3).^2)+2); %Función que representa la temperatura
surf(x,y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)
}}

== Gradiente de temperatura ==
Como se ha mencionado con anterioridad, el campo escalar que gobierna la temperatura en la placa es <math> T(x,y)=log((x-3)^2)+2)</math>.
 
El gradiente de un campo vectorial se define como: 
<center><math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial x} & \frac{\partial u_{1} }{\partial y} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial x} & \frac{\partial u_{2} }{\partial y} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial x} & \frac{\partial u_{3} }{\partial y} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> </center> 
Al no tener el campo componentes en y ni en z, la matriz resulta: 
<center> <math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial x} & 0 & 0\\
\frac{\partial u_{2} }{\partial x} & 0 & 0 \\
\frac{\partial u_{3} }{\partial x} & 0 & 0
\end{pmatrix} </math> </center>

De lo que resulta: <math>\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=>\frac{2(x-3)}{(x-3)^2+2} \vec i</math>
 
A continuación se muestra una representación del gradiente de la temperatura, como se puede apreciar este es ortogonal a las curvas de nivel de la temperatura ya que marca la dirección de máximo crecimiento de la misma.
Para este caso se observa que la temperatura crece en la dirección -<math>\vec i</math>
[[Archivo:Gradiente de temperatura_1.png]]
 
El código empleado para representar el gradiente ha sido el siguiente:
{{matlab|codigo=
h=0.2; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=pi/4:h:3*pi/4; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
x=ro.*cos(theta);
y=ro.*sin(theta);
z=y*0;;
T=log(((x-3).^2)+2); %Función de temperatura
grad=(2.*(x-3))./((x-3).^2+2); %Gradiente de temperatura
hold on
quiver(x,y,grad,0.*grad); %Se dibuja el gradiente
contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
axis([-3,3,-1,3]);
title("Gradiente de temperatura");
hold off
}}

== Campo de vectores ==
Se va a realizar una representación de el campo de desplazamientos que viene gobernado por la expresión: <math> \vec r(ρ,θ)= \frac{log(ρ)}{2}(sin(2θ-\frac{π}{2}))</math>
[[Archivo:Campo de desplazamientos sobre la curva.png]]
 
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
h=0.1; %Paso de muestreo
r=1:h:2; %Discretización del radio
t=1/4*pi:pi/50:3/4*pi; %Discretización del ángulo
[ro,theta]=meshgrid(r,t); %Generación de la retícula rectangular
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*(-sin(theta)); %Paso de la coordenada x del campo a cartesianas
FY=(log10(ro)./2).*(sin(2.*(theta)-(pi/2))).*cos(theta); %Paso de la coordenada y del campo a cartesianas
quiver(X,Y,FX,FY); %Dibuja el campo vectorial
axis equal
axis([-3,3,-1,3]);
view(2) %Vista en planta
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')
}}
 
== Comparación de la placa antes y después del desplazamiento ==
En este apartado se va a realizar una comparación visual de la placa antes y después de sufrir el desplazamiento.
[[Archivo:Placa antes y después.png]]
[[Archivo:Comparación de placas.png]] 
{{matlab|codigo=
El código empleado para llevar a cabo la comparación se expone a continuación:
h=0.2;%paso de muestreo
r=1:h:2;
t=pi/4:h:((3*pi/4)+0.2);
[ro,theta]=meshgrid(r,t); %creamos el mallado
x=ro.*cos(theta);
y=ro.*sin(theta);
%sólido antes del desplazamieneto
subplot(1,2,1);
i=mesh(x,y,x.*0);
view(2);
set(i,'EdgeColor','r')
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa no desplazada');
%sólido despues del desplazamiento
subplot(1,2,2);
%pasamos el campo de desplazamiento a cartesianas
A=-(log10(ro)./2).*(sin(2.*theta-(pi/2)).*(sin(theta)));
B=(log10(ro)./2).*(cos(theta)).*(sin(2.*theta-(pi/2)));
%Este campo de desplazamiento no desplaza prácticamente la placa, si no que concentra la masa de la misma en el centroide.
X=x+A;
Y=y+B;
j=mesh(X,Y,0.*X);
view(2);
set(j,'EdgeColor','b');
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');;
ylabel('Eje y');
title('Placa desplazada');
%Comparación dde ambas placas
hold on
figure
j=mesh(X,Y,X*0);
axis equal
axis([-3,3,-1,3]);
view(2);
set(j,'EdgeColor','b');
title('Comparación de ambas placas: Roja (inicial), azul (final)');
hold on
i=mesh(x,y,x.*0);
set(i,'EdgeColor','r')
hold off
}}
== Rotacional ==
El rotacional de un campo vectorial permite observar la tendencia que tiene dicho campo a inducir rotación alrededor de un punto de la placa.
El rotacional, calculado en coordenadas cilíndricas viene determinado por la siguiente expresión:
 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|</math>. 
Si se aplica al campo de desplazamientos resulta: 
<math>\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ·\frac{log(ρ}{2}·sin(2θ-\frac{π}{2}) & 0 \end{matrix}\right|=> -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}</math>. 
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
[[Archivo:Rotacional (2).png]] 
Como puede observarse, el rotacional es mayor en los puntos centrales de la placa, es por ello por lo que en la vista en 3 dimensiones la placa resulta doblada por esa zona. 
Del mismo modo, el programa nos indica que el valor máximo que alcanza el rotacional es 0.1504. 
El código empleado para dibujar el rotacional del campo de desplazamientos ha sido el siguiente:
{{matlab|codigo=
%Se crean los intervalos de las variables
h=50;
r=linspace(1,2,h);
t=linspace(pi/4,3*pi/4,h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=-(cos(2.*theta)./(2.*ro)).*log10((ro)+1)
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
shading interp
axis equal;
view(2);
colorbar;
title('Rotacional en 2D')
xlabel('Eje x');
ylabel('Eje y');
%Gráfica en 3D
subplot(1,2,2);
surf(x,y,rotacional);
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
Maximo=max(max(rotacional));
fprintf('el valor del rotacional máximo alcanzado es de: %1.4f \n',Maximo);
}}

==Tensor de deformaciones==
A continuación, se va a analizar el comportamiento de la placa al estar sometida a una serie de tensiones en un medio elástico lineal, isótropo y homogéneo. Se va a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento ondulatorio de los puntos de la misma dado por el siguiente campo vectorial: 
<math>\vec u(x,y,t)=\vec asin((k(\vec d \vec r -vt))</math>
 
Donde los desplazamientos permiten escribir el tensor de tensiones como: 
<math>σ=λ·\nabla·\vec u I+2με</math> 
Para ello, se va a definir 
<math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> 
Como la parte simétrica del tensor gradiente <math>\vec u</math> conocido como tensor de deformaciones. 
λ y μ son los llamados coeficientes de Lamé que por las características del material se conoce que son igual a 1.
Se define el tensor de deformaciones <math>\epsilon(\vec u(ρ,θ))</math> como la parte simétrica del tensor <math> \nabla{\vec u(ρ,θ)} </math>, llamado también gradiente del campo de desplazamientos <math>\vec u(ρ,θ) </math>:
: <math>\epsilon(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}</math>.
 
Lo primero será calcular <math>\nabla{\vec u(ρ,θ)}</math> que resulta: 
<math>∇\vec{u} =
\begin{pmatrix}
\frac{\partial u_{1} }{\partial ρ} & \frac{\partial u_{1} }{\partial θ} & \frac{\partial u_{1} }{\partial z} \\
\frac{\partial u_{2} }{\partial ρ} & \frac{\partial u_{2} }{\partial θ} & \frac{\partial u_{2} }{\partial z} \\
\frac{\partial u_{3} }{\partial ρ} & \frac{\partial u_{3} }{\partial θ} & \frac{\partial u_{3} }{\partial z}
\end{pmatrix} </math> 
 
 
::*<math>\frac{\partial \vec u}{\partial ρ}=(\frac{1}{2ρ})sen(2θ-\frac{π}{2}) \vec e_θ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial θ}= log(ρ)cos(2θ-\frac{π}{2}) \vec e_θ - (\frac{log(ρ)}{2})sen(2θ-\frac{π}{2}) \vec e_ρ</math>
 
 
::*<math>\frac{\partial \vec u}{\partial z}= 0 \vec e_z</math>
 
 
El gradiente del campo de deformaciones resulta: 
<math>\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2})& 0\\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
A continuación se calcula la traspuesta de la matriz del gradiente del campo de deformaciones: 
<math>(\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{1}{2ρ})sen(2θ-\frac{π}{2}) & 0 \\ -(\frac{log(ρ)}{2ρ})sen(2θ-\frac{π}{2}) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>.
 
 
Aplicando la fórmula <math>ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}</math> el tensor de deformaciones será igual a: 

<math>\epsilon(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{4ρ})(1 - log(ρ)) & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ 0 & 0 & 0 \end{matrix}\right)</math>

 
 

===Tensiones normales===
====Tensiones normales en la dirección del eje <math>\vec e_ρ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_ρ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math> 
 
Que resulta: <math>(\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e rho.png|500px||left|]] [[Archivo:3drho_1.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
a=(log10(ro).*sin(2.*theta))./(ro);%elemento(1,1,1);
figure(1);
surf(x,y,a);
shading interp
axis equal
title('Tensión normal en la dirección e_ρ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_θ</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_θ</math> es el siguiente:
 
<math> |(σ·\vec e_ρ))\vec e_ρ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right) ·\left ( \begin{matrix} 1 \\0\\0\end{matrix} \right ) ·\vec e_ρ</math>
 
 
Que resulta: <math> (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
 
<center>[[Archivo:Tension normal en e theta.png|500px|left]] [[Archivo:3drho.png|500px|right]]</center> 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

El código empleado ha sido el siguiente:
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
b=((log10(ro))./(ro))+2.*(log10(ro)).*sin(2.*theta);%elemento(2,2,2)
figure(1);
surf(x,y,b);
shading interp
axis equal
title('Tensión normal en la dirección e_θ')
view(2);
colorbar;
}}

====Tensiones normales en la dirección del eje <math>\vec e_z</math>====
El módulo de la tensión normal aplicada en la placa según la dirección del eje <math>\vec e_z</math> es el siguiente:
<math>|(σ·\vec e_θ))\vec e_θ | = \left ( \begin{matrix} (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & 0 \\ (\frac{sen(2θ-\frac{π}{2})}{2ρ})(1 - log(ρ)) & (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) & 0 \\ & \\0 & 0 & (\frac{log(ρ)}{ρ})cos(2θ-\frac{π}{2}) \end{matrix} \right ) ·\left ( \begin{matrix} 0 \\1\\0\end{matrix} \right ) ·\vec e_θ</math>
 
 
Que resulta: <math> (\frac{3log(ρ)}{ρ})cos(2θ-\frac{π}{2}) </math>
[[Archivo:Tension normal en e z.png|500px|left]] [[Archivo:3dz.png|500px|right]] 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
El código empleado ha siso el siguiente: 
{{matlab|codigo=
r=1:0.1:2;
t=pi/4:pi/50:3*pi/4;
[ro,theta]=meshgrid(r,t);
clf
x=ro.*cos(theta);
y=ro.*sin(theta);
z=zeros(size(x));
%se calculan las componentes de la matriz sigma
c=(log10(ro).*sin(2.*theta))./(ro);%elemento(3,3,3)
figure(1)
surf(x,y,c);
shading interp
axis equal
title('Tensión normal en la dirección e_k')
view(2);
colorbar;
}}

===Tensiones tangenciales respeto al plano ortogonal a <math>\vec e_ρ</math>===
===Tensión de Von Mises===
La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. 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 teoría expone que un material dúctil comenzará a ceder en una ubicación cuando la tensión de von Mises sea igual al límite de tensión. En resumen 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. 
Esta viene determinada por 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> 
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).
Para este supuesto, la tensión de Von Mises representará la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
 
A continuación se muestra la representación de dicha tensión:
[[Archivo:Von mises22b.png|centre|]]
 
El código empleado para la misma es:

[[Categoría:Teoría de Campos]]

Archivo:Von mises22b.png

2022-12-08T08:08:40Z

Sergio puente: