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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo 22-B
Asignatura	Teoría de Campos
Curso	2022-23
Autores	Alejandro Seises López Sergio Puente Puente Abel Larrosa Morales Víctor Sepúlveda Fernández
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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

2 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.
El código empleado para dibujar la placa es el siguiente:

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

3 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.
El código empleado para dibujar las curvas de nivel del campo de temperaturas es el siguiente:

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

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)

4 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}=\gt\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]
El código empleado para representar el gradiente ha sido el siguiente:

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

5 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]
El código empleado para dibujar el campo de desplazamientos ha sido el siguiente:

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

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

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

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:

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

8 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|=\gt -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}[/math].
A continuación se muestra una tanto vista del rotacional visto en planta como en 3 dimensiones.
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:

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

9 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]

10 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]

10.1 Tensiones normales

10.1.1 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]

El código empleado ha sido el siguiente:

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;

10.1.2 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]

El código empleado ha sido el siguiente:

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;

10.1.3 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]
















El código empleado ha siso el siguiente:

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;

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

El código empleado ha sido el siguiente:

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

10.3 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:
El código empleado ha sido el siguiente:

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

11 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é equivalente:

[math] \vec F=\frac{\partial^2 \vec u}{\partial t^2}-\mu \Delta \vec u-(\lambda+\mu)\nabla(\nabla \cdot \vec u) \, [/math]

Por las propiedades del material se conoce que μ y λ=1, por lo que la ecuación resulta:

[math] \vec F=-\Delta \vec u-2\nabla(\nabla \cdot \vec u) [/math]

A continuación se muestra una representación gráfica del campo de fuerzas:

El código empleado ha sido el siguiente:

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