Estudio de las transformaciones de un cuarto de anillo (Grupo 13B)

Trabajo realizado por estudiantes
Título	Deformaciones en un cuarto de anillo.
Asignatura	Teoría de Campos
Curso	2022-23
Autores	Felipe Delgado Rubio, Carlos Reillo González, Jorge Granadino Aranda, Diego de la Flor Bravo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa plana definida por un anillo centrado en el origen comprendido entre los radios 1 y 2, y por el plano \(y≥|x|\). Para realizar el estudio dispondremos de la función de temperatura, dada en coordenadas cartesianas:
[math]T(x,y)=log((x-3)^2+2)[/math].

A su vez, también dispondremos de la función de campo, dada en coordenadas cilíndricas:
[math]\vec u(ρ,θ)=\frac{log(ρ)}{2}sin(2θ-\frac{π}{2}) ~[/math]\(\vec e_θ\).

1 Mallado

Nuestro anillo esta regido por varias condiciones; el radio debe estar contenido entre 1 y 2, y el anillo está comprendido en el plano \(y≥|x|\). Estas condiciones en coordenadas cilíndricas se representan como; (ρ,θ) ∈ [1,2]×[[math]\frac{π}{4}[/math],[math]\frac{3π}{4}[/math]]. Tomando como ejes para la representación del mallado del anillo, \((x, y) ∈ [-3; 3] × [−1; 3]\) y tomando el paso de muestreo h=[math]\frac{2}{10}[/math] en los ejes x e y, permitirá realizar una representación exacta de la placa en el espacio 2-D.

Representación del mallado

h=2/10;  %Paso de muestreo
rho=1:h:2;  %Asignación de las variables rho y theta 
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
mesh(x,y,0.*x)
view(2)
axis([-3,3,-1,3])  %Asignación del eje
xlabel('Eje X')
ylabel('Eje Y')
title('Representación placa en 2D')

2 Curvas de nivel de la temperatura

Las curvas de nivel son la representación de la temperatura sobre la placa, sabiendo que la función temperatura es;

[math]T(x,y)=ln((x−3)^2 +2)[/math]

Representación de las curvas de nivel

h=2/10;  %Paso de muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
T=log((x-3).^2+2);  %Función de temperaturas
hold on
subplot(1,2,1)
surf(x,y,T)  %Representación el campo escalar de temperaturas
view(2)
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
colorbar
subplot(1,2,2)
contour(x,y,T,20/h)
title('Curvas de Nivel')
colorbar
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
hold off

La temperatura máxima es 3,0674ºC,

Maximo=max(max(T));  %Asigna el numero máximo de la función T a la variable 'Maximo'

3 Gradiente de la temperatura

El gradiente de la temperatura representa la dirección en la que varía más rápido. El gradiente se obtiene mediante la siguiente fórmula; [math]\nabla T(x,y)=\frac{∂T}{∂x}\vec i +\frac{∂T}{∂y}\vec j[/math].

Lo que hace que el gradiente de [math]T(x, y)=ln((x−3)^2 +2)[/math], sea:

[math]\nabla T(x,y)=\frac{2(x-3)}{(x-3)^2+2} \vec i[/math]

Para una mejor interpretación, el gradiente será representado en 2D y en 3D

Representaciones del gradiente de la temperatura

h=2/10;  %Paso de Muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]= meshgrid(rho,tt);  %Creación del mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
G=(2.*(x-3)./((x-3).^2)+2);  %Gradiente de la temperatura
T=log((x-3).^2 + 2);  %Función de temperatura
subplot(1,2,2)  %Gráfica en 3D
hold on
surf(x,y,T)
quiver3(x,y,T,G,0*y,0*T)
axis vis3d
view(3)
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Gradiente 3D')
hold off
subplot(1,2,1)  %Gráfica en 2D
view(2)
hold on
contour(x,y,T,60)
ti=inline('(2.*(x-3))./(((x-3).^2)+2)','x','y');
tj=inline('0.*x','x','y');
Tx=ti(x,y);
Ty=tj(x,y);
quiver(x,y,Tx,Ty)
axis([-3,3,-1,3])
title('Gradiente de la temperatura en 2D')
hold off

4 Desplazamiento producido por [math]\vec u~[/math]

El desplazamiento que se produce en la placa es el que ejerce [math]\vec u(ρ,θ)[/math] sobre la placa.

h=2/10;  %Paso de muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
subplot(1,2,1)  %Placa en el inicio
m=mesh(x,y,0*x);
view(2)
set(m,'EdgeColor','r')  %Color de la placa de inicio
axis([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Placa no desplazada')
subplot(1,2,2)  %Placa desplazada
A=(sin(TT).*cos(2.*TT).*log(RR))./2;
B=(-cos(TT).*cos(2.*TT).*log(RR))./2;
X=x+A;
Y=y+B;
n=mesh(X,Y,0*X);
view(2)
set(n,'EdgeColor','b')  %Color de la placa desplazada
axis([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Placa desplazada')
figure  %Comparación de ambas placas
n=mesh(X,Y,0*X);
axis([-3,3,-1,3])
view(2)
set(n,'EdgeColor','b')  %Color de la placa desplazada
title('Desplazamiento de la placa')
hold on
m=mesh(x,y,0*x);
set(m,'EdgeColor','r')  %Color de la placa de inicio
hold off

La placa antes y después de la acción del campo

Comparación entre antes (rojo) y después del desplazamiento (azul)

Se puede observar que la superficie de la placa no cambia de forma significativa; sin embargo, se puede apreciar como actúa el campo sobre ella.

5 Divergencia

La divergencia es un operador vectorial que opera sobre un campo vectorial, produciendo un campo escalar. Representa los cambios de volumen local originados por el vector velocidad, en nuestro caso el vector velocidad es; [math]\vec u(ρ,θ)=\frac{ln(ρ)}{2}sin(2θ-\frac{π}{2})\vec e_θ [/math].

La divergencia en cilíndricas se calcula;

[math]\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho·u_ρ})+\frac{\partial}{\partial{\theta}}(u_θ))] [/math]

Por lo tanto, la divergencia de nuestro campo vectorial es;

[math]\nabla\cdot\vec u=\frac{ln(ρ)}{ρ}sin(2θ)[/math]

Para entender de mejor manera la divergencia, se graficara en 2D y 3D

Representaciones de la divergencia del campo vectorial

h=2/10; %Paso de Muestreo
rho=1:h:2;  %Asignación de las variables rho y theta 
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt); %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
DIVu=(1./RR).*(log(RR).*sin(2.*(TT))); %Divergencia de u
subplot(1,2,1)  %Gráfica 2D
surf(x,y,DIVu)
shading interp
view(2)
axis([-3,3,-1,3])
colorbar
title('Divergencia 2D')
xlabel('Eje X')
ylabel('Eje Y')
subplot(1,2,2)  %Gráfica 3D
view(2)
surf(x,y,DIVu)
shading interp
axis([-3,3,-1,3])
colorbar
axis vis3d
title('Divergencia 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')

La divergencia es máxima en el punto ([math]\sqrt{2}[/math],[math]\sqrt{2}[/math]); donde vale 0.3466, mínima en el punto ([math]-\sqrt{2}[/math],[math]\sqrt{2}[/math]); donde vale -0.3466, y nula en (0,0). Estos puntos se pueden apreciar en la gráfica.

6 Rotacional

El rotacional es un operador vectorial que actúa sobre campos vectoriales, muestra la tendencia de un campo vectorial a producir rotación en un punto. La expresión del rotacional calculada en coordenadas cilíndricas se define por la siguiente expresión;

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

Aplicándola sobre el campo conocido, [math]\vec u(ρ,θ) [/math], devuelve la expresión;

[math]\nabla×\vec u(ρ,θ)=-\frac{cos(2θ)(ln(ρ)+1)}{2ρ} \vec e_z[/math]

Rotacional

h=2/10;  %Paso de Muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
ROTu=(-cos(2.*TT).*(log(RR)+1))./(2.*RR) ;  %Rotacional de u
subplot(1,2,1)  %Representación en 2D
surf(x,y,ROTu)
shading interp
axis([-3,3,-1,3])
view(2)
title('Rotacional en 2D') 
xlabel('Eje X')
ylabel('Eje Y')
subplot(1,2,2)  %Representación en 3D
surf(x,y,ROTu)
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')

El punto que sufre mayor rotacional es el punto (0,1).

7 Tensor de tensiones

El objeto del apartado es determinar σ. Una vez definida la parte simétrica del tensor gradiente de u conocido como tensor de deformaciones. siendo este ε(u). Entonces conoceremos el tensor de tensiones definido como σ. Sabemos que σ se puede obtener a partir de la fórmula:

centro

Tomamos los coeficientes de Lamé (λ y μ) ambos como iguales a 1, y el tensor de deformaciones como la suma del gradiente de \(\vec u\) con su matriz traspuesta divididos por 2. EL resultado de sigma lo podemos ver incluido en el código del programa.

Después podremos representar las tensiones normales en la dirección que marca el eje i, es decir [math]\vec i [/math] = \(\vec g_{i}\) * σ * \(\vec g_{i}\), las tensiones normales en la dirección que marca el eje j, es decir [math]\vec j [/math] = \(\vec g_{j}\) * σ * \(\vec g_{j}\) y las correspondientes al eje k, es decir [math]\vec k[/math] = \(\vec g_{k}\) * σ * \(\vec g_{k}\). Resultan las siguientes gráficas:

h=2/10;  %Paso de Muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
a=log(RR).*sin(2.*TT)./RR;  %Elemento (1,1) y (3,3) de la matriz
b=(1./2.*RR).*(RR.*log(RR).*cos(2.*TT)-cos(2.*TT)); %Elemento (1,2) y (2,1) de la matriz
c=2.*log(RR).*sin(2.*TT)./RR;  
d=a+c  %Elemento (2,2) de la matriz
subplot(1,3,1) %Tensión normal en e_rho
surf(x,y,a)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e sub rho')
subplot(1,3,2)  %Tensión normal en e_theta
surf(x,y,d)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e sub theta')
subplot(1,3,3)  %Tensión normal en e_z
surf(x,y,a)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e_z')
figure
subplot(1,3,1)  %Tensión normal en e_rho
surf(x,y,a)
axis equal
shading interp
colorbar
title('Tensión normal eje e sub rho')
subplot(1,3,2)  %Tensión normal en e_theta
surf(x,y,d)
axis equal
shading interp
colorbar
title('Tensión normal eje e sub theta')
subplot(1,3,3)  %Tensión normal en e_z
surf(x,y,a)
axis equal
shading interp
colorbar
title('Tensión normal eje e_z')

Representación de las tensiones normales a los ejes

Vista plano XY

8 Tensiones tangenciales

Habiendo calculado σ en el apartado anterior, se determinan las tensiones tangenciales respecto al plano ortogonal a \(\vec e_ρ\). Únicamente quedan representadas aquellas tensiones que no son nulas. Se presenta el módulo de la tensión tangencial como:

[math]|σ·\vec e_ρ-(\vec e_ρ·(σ·\vec e_ρ))\vec e_ρ|[/math]

Representaciones de las tensiones tangenciales

h=2/10;  %Paso de muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
tang=((-cos(2.*TT)./2.*RR).*(1-log(RR)));  %Función de la tensión tangencial
subplot(2,1,1)  %Representación en 2D
surf(x,y,tang)
shading interp
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(2)
colorbar
title('Tensión tangencial a la dirección del eje e sub rho')
subplot(2,1,2)  %Representación en 3D
surf(x,y,tang)
shading interp
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
colorbar
title('Tensión tangencial a la dirección del eje e sub rho')

9 Tensión de Von Mises

Para este apartado tratamos la Tensión de Von Mises. Se presenta como una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico. Emplearemos la siguiente fórmula, que toma como valores de σ(1), σ(2) y σ(3) los autovalores de la matriz σ (también conocidos como tensiones principales). Para completar el apartado se pinta la Tensión de Von Mises y se señala el punto en el que alcanza el mayor valor

Representación de la tensión de Von Mises

h=2/10;  %Paso de muestreo
rho=1:h:2;  %Asignación de las variables rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Malldo
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
sig= [];  %Creación de matriz sigma con tantos columnas como la longitud de theta y tantas columnas como la longitud de rho
VM=zeros(length(tt),length(rho));
a = @(R,T) (log(R)./R).*sin(2.*T);  %Fórmulas de las componentes de sigma
b = @(R,T) (-cos(2.*T)./(2.*R)).*(1-log(R));
c = @(R,T) (3.*log(R)./R).*sin(2.*T);

for i=1:length(tt)  %Definición de las componentes de la matriz de Von Mises
  for j=1:length(rho)
    sig(1,1)= a(RR(i,j),TT(i,j));
    sig(1,2)= b(RR(i,j),TT(i,j));
    sig(1,3)= 0;
    sig(2,1)= b(RR(i,j),TT(i,j));
    sig(2,2)= c(RR(i,j),TT(i,j));
    sig(2,3)= 0;
    sig(3,1)= 0;
    sig(3,2)= 0;
    sig(3,3)= a(RR(i,j),TT(i,j));
    
    [v,u]=eig(sig);  %Cálculo de autovalores
    VM(i,j)=sqrt(((u(1,1)-u(2,2))^2+(u(2,2)-u(3,3))^2+(u(3,3)-u(1,1))^2)/2);  %Aplicación de la fórmula
  end
end
subplot(1,3,1)  %Representación en 2D en el eje XOY
surf(x,y,VM)
shading interp
view(2)
xlabel('Eje x')
ylabel('Eje y')
axis equal
colorbar;
title('XOY')
subplot(1,3,2)  %Representación en 3D
surf(x,y,VM)
shading interp
axis equal
axis vis3d
colorbar;
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
title('Tensión de Von Mises')
subplot(1,3,3)  %Representación en 2D en el eje XOZ
hold on
m=VM(1:size(VM,1),1);
t=x(1:length(x),3);
surf(x,y,VM)
view(0,0)
shading interp
MAX = max(m);
for k = 1:length(m)  %Cálculo del punto máximo
  if m(k) == MAX
    plot3(t(k),0,MAX,'xr','markersize',10)  %Representación del punto máximo
  endif
endfor
axis equal
colorbar
title('XOZ')
xlabel('Eje x')
zlabel('Eje z')
txt = ['Máxima tensión:' num2str(MAX)];  %Texto insertado sobre el dibujo indicando el máximo
text(-0.35,2,1.25,txt)
hold off

10 Campo de fuerzas [math] \vec F [/math].

Las fuerzas que actuan sobre la placa y que causan el desplazamiento observado, se calculan con la ecuacion de la elasticidad lineal [math]\vec F [/math]:
[math] \vec{F}=\frac{\partial^2\vec{u} }{\partial t^2}-\mu \Delta \vec{u} -(\lambda +\mu )\nabla(\nabla\cdot\vec{u})[/math]
La ecuación añadida nos da como resultado un campo vectorial que representa unas fuerzas aplicadas sobre la placa. Para cada punto del campo vectorial, [math]\vec F [/math] devuelve un vector que representa la fuerza aplicada en dicho punto.
El Laplaciano de un campo vectorial, operación que aparece en la ecuación, se calcula mediante la siguiente fórmula:
[math]\Delta \vec{u}=\Delta (u_1\vec{i}+u_2\vec{j}+u_3\vec{k})=\Delta u_1\vec{i}+\Delta u_2\vec{j}+\Delta u_3\vec{k}[/math]
También, tomamos los coeficientes de Lamé (λ y μ) ambos iguales a 1.

Representación del campo de fuerzas

h=2/10;  %Paso de muestreo
rho=1:h:2;  %Asignación de los valores de rho y theta
tt=linspace(pi/4,3*pi/4,pi/h);
[RR,TT]=meshgrid(rho,tt);  %Mallado
x=RR.*cos(TT);  %Cambio a coordenadas cilíndricas
y=RR.*sin(TT);
i=((2.*(log(RR)-1).*sin(2.*TT))./(RR.^2))-(log(RR).*(sin(TT)-9.*sin(3.*TT))./((2.*RR).^2));  %Fórmulas de la Fuerza en i, j
j=-((4.*cos(2.*TT).*log(RR))./(RR.^2))-(log(RR).*(cos(TT)+9.*sin(3.*TT))./((2.*RR).^2));
hold on
mesh(x,y,0.*x)
quiver(x,y,i,j)
axis ([-2,2,-0,3])
title('Campo de fuerzas')
hold off