Análisis y visualización de transformaciones en un cuarto de anillo

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Análisis y visualización de transformaciones en un cuarto de anillo (Grupo 18)
Asignatura Teoría de Campos
Curso 2023-24
Autores Sara López Caro
Luis Melgar Beltrán
Pablo Sancho González
Mª Eugenia Sanz de Ojeda
Clara Sanz Riomoros
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


Se considera una placa plana definida por un cuarto de un anillo centrado en el origen y comprendido entre los radios 1 y 2, y por el plano \(y≥|x|/2\). El grupo hará un análisis de diferentes transformaciones y aplicaciones que se le pueden ejercer a dicha placa.

Para la realización del estudio, se dispone de la función de temperatura, dada en coordenadas cartesianas:
[math]T(x,y)=\cos((y-3)^2+x)[/math].

Y de la función de campo, dada en coordenadas cilíndricas:
[math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1}sin(2θ-\frac{π}{2}) ~[/math]\(\vec e_θ\).

1 Mallado del sólido

En primer lugar, se visualizará la placa considerada. Con las características descritas en el enunciado:

  • Una placa plana que ocupa un cuarto de un anillo circular centrado en el origen y comprendido entre los radios 1 y 2, en el plano y ≥ |x|/2.
  • Dibujado con los ejes en el cuadrado [−3, 3] × [−1, 3]
h=2/10;  %Paso de muestreo
ro=1:h:2;  %Asignación de las variables rho y theta
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]=meshgrid(ro,zeta);  %Mallado
x=RR.*cos(Z);  %Cambio a coordenadas cilíndricas
y=RR.*sin(Z);
mesh(x,y,0.*x)
view(2)
axis([-3,3,-1,3])  %Asignación del eje
xlabel('Eje X')
ylabel('Eje Y')
title('Mallado de la placa')

Mallado de la placa


2 Curvas de nivel de la temperatura

El enunciado ofrece una función de temperatura, en la que los valores de dicha temperatura varían según el punto que se estudie. La función de temperatura dada es, en coordenadas cartesianas:
[math]T(x,y)=cos((y-3)^2+x)[/math].

2.1 Visualización curvas de nivel

Se dibujarán sobre la placa las distintas curvas de nivel que representan los lugares de la placa que comparten temperatura. Para esto se utilizará el comando de Matlab: Contour. 

h = 2/10;  % Paso de muestreo
ro = 1:h:2;  % Asignación de las variables ro y zeta
zeta = linspace(atan(1/2), pi-atan(1/2), pi/h);
[RR, Z] = meshgrid(ro, zeta);  % Mallado
x = RR.*cos(Z);  % Cambio a coordenadas cilíndricas
y = RR.*sin(Z);
T = cos((y-3).^2 + x);  % Función de temperaturas T(x, y) = cos((y − 3)2 + x)
[max_temp, i] = max(T(:)); % Máximo de la función de temperatura
[coord_fila, coord_col] = ind2sub(size(T), i);
coord_max_temp = [x(coord_fila, coord_col), y(coord_fila, coord_col)];
% Graficación del resultado
hold on		% Superposición de gráficas
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,15/h)
title('Curvas de Nivel')
colorbar
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
hold off
disp(coord_max_temp)
disp(max_temp)
colorbar;


Matlab ofrece, con este código, las coordenadas del punto donde la temperatura es máxima y su temperatura. En este caso la temperatura máxima es 0.9999 ºC y se alcanza en (0.9754,0.6990) Representación del campo escalar de temperaturas y sus curvas de nivel

2.2 Visualización del gradiente de temperatura

El gradiente de una función de temperatura representa la tasa de cambio espacial de la temperatura en un punto específico, indicando la dirección y magnitud en la que la temperatura aumenta con mayor rapidez. Se visualizarán los vectores que describen esto. El gradiente de la temperatura se obtiene aplicando la fórmula :[math]\nabla T(x,y)=\frac{∂T}{∂x}\vec i +\frac{∂T}{∂y}\vec j[/math] en la función de temperatura dada:

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

Obteniéndose así la siguiente función:
[math] ∇T=-sin((y-3)^2+x)\vec i-(2y-6)sin((y-3)^2+x)\vec j [/math] 

h=2/10;  %Paso de Muestreo
ro=1:h:2;  %Asignación de las variables rho y theta
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]= meshgrid(ro,zeta);  %Creación del mallado
x=RR.*cos(Z);  %Cambio a coordenadas cilíndricas
y=RR.*sin(Z);
Gi=-sin(y.^2-6.*y+9+x);
Gj=-(2.*y-6).*sin(y.^2-6.*y+9+x);  %Gradiente de la temperatura
T = cos((y-3).^2 + x);  %Función de temperatura
subplot(1,2,2)  %Graficación en 3D
hold on
surf(x,y,T)
quiver3(x,y,T,Gi,Gj,0*T) 
view([75 25])
axis equal
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Gradiente 3D')
hold off
subplot(1,2,1)  %Graficación en 2D
view(2)
hold on
contour(x,y,T,60)
ti = @(x, y) -sin(y.^2-6.*y+9+x);
tj = @(x, y) -(2.*y-6).*sin(y.^2-6.*y+9+x);
Tx=ti(x,y);
Ty=tj(x,y);
quiver(x,y,Tx,Ty)
axis equal
axis([-3,3,-1,3])
title('Gradiente de la temperatura en 2D')
hold off


Representación del gradiente de temperatura sobre la placa en dos y tres dimensiones.

3 Energía Calorífica y Ley de Fourier

Según la Ley de Fourier, la energía calorífica viaja mediante la siguiente fórmula:

[math]Q=-k∇T[/math]
Donde k corresponde a la constante de conductividad térmica, que en este caso será igual a 1, y ∇T es el gradiente de la temperatura. En términos sencillos, la fórmula indica que la cantidad de calor transferido es proporcional a la constante de conductividad térmica y a la variación de temperatura. 

h=2/10;  %Paso de Muestreo
ro=1:h:2;  % Asignación de las variables ro y zeta
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]= meshgrid(ro,zeta);  %Creación del mallado
x=RR.*cos(Z);  % Cambio a coordenadas cilíndricas
y=RR.*sin(Z);
Gi=-sin(y.^2-6.*y+9+x);
Gj=-(2.*y-6).*sin(y.^2-6.*y+9+x);   % Gradiente de la temperatura
T=cos((y-3).^2 + x);  % Función de temperatura
hold on
mesh(x,y,0.*x)
ti= @(x,y) sin(y.^2-6.*y+9+x);
tj= @(x,y) (2.*y-6).*sin(y.^2-6.*y+9+x);
Qx=ti(x,y);
Qy=tj(x,y);
quiver(x,y,Qx,Qy)
axis equal
axis([-3,3,-1,3])
hold off
title('Movimiento Energía Calorífica')


Representación del movimiento de la energía calorífica.

4 Campo de vectores

Para realizar el campo de vectores en los puntos de mallado del sólido, se dispone de la función de campo:
[math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1}sin(2θ-\frac{π}{2}) ~[/math]\(\vec e_θ\). 

h=2/10;
ro=1:h:2;
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]=meshgrid(ro,zeta);
x=RR.*cos(Z);
y=RR.*sin(Z);
mesh(x,y,y*0); % Mallado del sólido
title('CAMPO DE VECTORES')
xlabel ('EJE X')
ylabel ('EJE Y')
view(2)
hold on
mx=(1/2).*exp(RR-1).*sin(2.*Z-(pi/2)).*cos(Z); %Campo de vectores
my=(1/2).*exp(RR-1).*sin(2.*Z-(pi/2)).*sin(Z);
quiver(x,y,mx,my);
axis equal
axis([-3,3,-1,3])
hold off


Ej4(AyV).png

5 Desplazamiento dado el campo de vectores

El desplazamiento que se produce en la placa es el que ejerce [math]\vec u(ρ,θ)[/math] sobre la placa. Por tanto, a las coordenadas que se tenían del cuarto de anillo, se le han sumado el desplazamiento de cada una de ellas, obteniendo unas nuevas coordenadas para cada uno de los puntos. Se compara en este apartado la disposición de la placa antes y después de dicho movimiento. 

h=2/10;
ro=1:h:2;
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]=meshgrid(ro,zeta);
x=RR.*cos(Z);
y=RR.*sin(Z);
% Principio
subplot(2,2,1)
surf(x,y,y*0);
axis([-3,3,-1,3]);
view(2)
title('SITUACIÓN INICIAL')
% final
subplot(2,2,2)
mx=(1/2).*exp(RR-1).*sin(2.*Z-(pi/2)).*cos(Z); %Cambio a cartesianas
my=(1/2).*exp(RR-1).*sin(2.*Z-(pi/2)).*sin(Z);
rx=x+mx;
ry=y+my;
surf(rx,ry,ry*0);
axis([-3,3,-1,4]); %Se agranda la cuadrícula hacia arriba para que se pueda percibir la transformación
view(2)
title('SITUACIÓN FINAL')
% comparación
subplot(2,2,3)
plot3(x,y,y*0);
hold on
plot3(rx,ry,ry*0);
axis([-3,3,-1,4]);
view(2)
title('COMPARACIÓN')
hold off


Ej5(AyV).jpg

6 Divergencia

La divergencia es un operador vectorial que opera sobre un campo vectorial, produciendo un campo escalar. Representa la cantidad por la cual el campo "se aleja" o "converge" en un punto dado. Con el campo:

[math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1}sin(2θ-\frac{π}{2}) ~[/math]\(\vec e_θ\).

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 tanto, en el caso del campo dado, la divergencia es:

[math] \nabla\cdot\vec u= \frac{e^{ρ-1}}{ρ}cos(2θ- \frac{π}{2}) [/math]

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

h=2/10;
ro=1:h:2;
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]=meshgrid(ro,zeta);
x=RR.*cos(Z);
y=RR.*sin(Z);
DIVu=(1./RR).*exp(RR-1).*cos(2*Z-pi/2); %Divergencia de u
subplot(1,2,1)  %Graficación 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)  %Graficación 3D
surf(x,y,DIVu)
shading interp
axis([-3,3,-1,3])
view(195,10)
colorbar
title('Divergencia 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')

% Encontrar el máximo, mínimo y nulo de la divergencia
[max_div, i] = max(DIVu(:));
[coord_fila, coord_col] = ind2sub(size(DIVu), i);
coord_max_div = [x(coord_fila, coord_col), y(coord_fila, coord_col)];
[min_div, i] = min(DIVu(:));
[coord_fila, coord_col] = ind2sub(size(DIVu), i);
coord_min_div = [x(coord_fila, coord_col), y(coord_fila, coord_col)];
[nulo_div, i] = min(abs(DIVu(:)));
[coord_fila, coord_col] = ind2sub(size(DIVu), i);
coord_nulo_div = [x(coord_fila, coord_col), y(coord_fila, coord_col)];

disp('La máxima divergencia es :')
disp(max_div)
disp(coord_max_div)
disp('La miníma divergencia es: ')
disp(min_div)
disp(coord_min_div)
disp(coord_nulo_div)


Ej6(AyV).jpg

La divergencia es máxima en el punto (1.4219,1.4065) con un valor de 1.3591.
La divergencia mínima se encuentra en (-1.4219,1.4065) con un valor de -1.3591.
Por último, la divergencia nula está en el punto ( 0.0000,1.0000).

En la gráfica en 3D se observa cómo cada punto cambia de ubicación y, por tanto, algunos de los diferenciales de área varían. Además, los puntos calculados anteriormente se pueden estudiar en la gráfica.

7 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} \\ \ u_ρ & \ ρu_θ & \ u_z \end{matrix}\right|[/math].

Por tanto, el rotacional del campo dado es:

[math] \nabla×\vec u = (1+ \frac{1}{ρ})\frac{e^{ρ-1}}{2}sen(2θ-\frac{π}{2}) \vec e_z [/math] 

h=2/10; 
ro=1:h:2;
zeta=linspace(atan(1/2),pi-atan(1/2),pi/h);
[RR,Z]=meshgrid(ro,zeta);
x=RR.*cos(Z);
y=RR.*sin(Z);
Rot_z=(1+1./RR)*(1/2).*exp(RR-1).*sin(2*Z-pi/2);
subplot(1,2,1)  %Gráfica 2D
surf(x,y,Rot_z)
shading interp
view(2)
axis equal
axis([-3,3,-1,3])
colorbar
title('Rotacional 2D')
xlabel('Eje X')
ylabel('Eje Y')
subplot(1,2,2)  %Gráfica 3D
surf(x,y,Rot_z)
shading interp
axis([-3,3,-1,3])
view(195,10)
title('Rotacional 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')

% Encontrar el máximo, mínimo y nulo del rotacional
[max_rot, i] = max(Rot_z(:));
[coord_fila, coord_col] = ind2sub(size(Rot_z), i);
coord_max_rot = [x(coord_fila, coord_col), y(coord_fila, coord_col)];

[min_rot, i] = min(Rot_z(:));
[coord_fila, coord_col] = ind2sub(size(Rot_z), i);
coord_min_rot = [x(coord_fila, coord_col), y(coord_fila, coord_col)];


disp('Máximo rotacional: ')
disp(max_rot)
disp(coord_max_rot)

disp('Mínimo rotacional: ')
disp(min_rot)
disp(coord_min_rot)


Izquierda Derecha

El punto con máximo rotacional es (0.0000,2.0000) con un valor de 2.7183.
El punto con mínimo rotacional es (1.7889,0.8944) con un valor de -1.6310.

8 Tensor de deformaciones

Si definimos el tensor deformaciones como:

[math] ε(ū)=(∇ū+∇ūt)/2 [/math].

Entendemos el tensor de tensiones σij:

[math] σ = λ∇·ū 1 + 2µε [/math].

λ, µ son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material, tomaremos λ = µ = 1.
Por tanto: [math] σ = λ∇·ū + ∇ū + ∇ūt [/math].

∇ū =

  • [math] ∇ū*\vec e_ρ = e^{ρ-1}*sen(2θ-\frac{π}{2})\vec e_θ [/math]
  • [math] ∇ū*\vec e_θ =\frac{e^{ρ-1}}{ρ}*cos(2θ-\frac{π}{2})\vec e_θ+\frac{e^{ρ-1}}{2}*sen(2θ-\frac{π}{2})(-\vec e_ρ)[/math]
  • [math] ∇ū*\vec e_z = 0*\vec e_θ + 0*\frac{∂\vec e_θ}{∂\vec e_z} = 0 [/math]


[math] ∇ū = \begin{pmatrix} \ 0 & \frac{-1}{2}*e^{ρ-1}*sen(2θ-\frac{π}{2}) & \ 0 \\ \ e^{ρ-1}*sen(2θ-\frac{π}{2}) & \frac{e^{ρ-1}}{ρ}*cos(2θ-\frac{π}{2}) & \ 0 \\ \ 0 & \ 0 & \ 0 \end{pmatrix}[/math]


[math] ∇ūt = \begin{pmatrix} \ 0 & \ e^{ρ-1}*sen(2θ-\frac{π}{2}) & \ 0 \\ \frac{-1}{2}*e^{ρ-1}*sen(2θ-\frac{π}{2}) & \frac{e^{ρ-1}}{ρ}*cos(2θ-\frac{π}{2}) & \ 0 \\ \ 0 & \ 0 & \ 0 \end{pmatrix}[/math]


[math] σ = \begin{pmatrix} \frac{e^{ρ-1}}{ρ}*cos(2θ-\frac{π}{2}) & \frac{e^{ρ-1}}{2}*sen(2θ-\frac{π}{2}) & \ 0 \\ \frac{e^{ρ-1}}{2}*sen(2θ-\frac{π}{2}) & \frac{3*e^{ρ-1}}{ρ}*cos(2θ-\frac{π}{2}) & \ 0 \\ \ 0 &\ 0 & \frac{e^{ρ-1}}{ρ}*cos(2θ-\frac{π}{2}) \end{pmatrix}[/math] 

h = 2/10;  % Paso de muestreo
ro = 1:h:2;  % Asignación de las variables rho y theta
zeta = linspace(atan(1/2), pi-atan(1/2), pi/h);
[RR, Z] = meshgrid(ro, zeta);  % Mallado
x = RR.*cos(Z);  % Cambio a coordenadas cilíndricas
y = RR.*sin(Z);
a=(1./RR).*exp(RR-1).*cos(2.*Z-(pi/2)); %Elemento (1,1) y (3,3) de la matriz
b=(1/2).*exp(RR-1).*sin(2.*Z-(pi/2)); %Elemento (1,2) y (2,1) de la matriz
c=3.*a; %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_r_o')
subplot(1,3,2)  %Tensión normal en e_theta
surf(x,y,c)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e_z_e_t_a')
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_r_o')
subplot(1,3,2)  %Tensión normal en e_theta
surf(x,y,c)
axis equal
shading interp
colorbar
title('Tensión normal eje e_z_e_t_a')
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')


Ej8.1 page-0001.jpg Ej8.2 page-0001.jpg


9 Tensiones tangenciales

Para calcular las tensiones tangenciales respecto a el plano otrogonal a i se aplica la siguiente fórmula:

[math] |σ ·~i − (~i · σ ·~i)~i| [/math].

Por tanto, en el caso del tensor de tensiones con el que se está trabajando:

[math] |σ ·~i − (~i · σ ·~i)~i| = \frac{e^{ρ-1}}{2}*sen(2θ-\frac{π}{2}) [/math]. 

h = 2/10;  % Paso de muestreo
ro = 1:h:2;  % Asignación de las variables rho y theta
zeta = linspace(atan(1/2), pi-atan(1/2), pi/h);
[RR, Z] = meshgrid(ro, zeta);  % Mallado
x = RR.*cos(Z);  % Cambio a coordenadas cilíndricas
y = RR.*sin(Z);
ti=(1/2)*exp(RR-1).*sin(2*Z-(pi/2));  %Función de la tensión tangencial
subplot(2,1,1)  %Representación en 2D
surf(x,y,ti)
shading interp
axis equal
axis([-3,3,-1,3])
xlabel('x')
ylabel('y')
zlabel('z')
view(2)
colorbar
title('Tensión tangencial respecto al plano ortogonal a i')
subplot(2,1,2)  %Representación en 3D
surf(x,y,ti)
shading interp
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
colorbar
title('Tensión tangencial respecto al plano ortogonal a i')


Representaciones de las tensiones tangenciales Representaciones de las tensiones tangenciales


10 La tensión de Von Mises

La tensión de Von Mises se define por la fórmula:

[math] σV M = sqrt (((σ1 − σ2)^2 + (σ2 − σ3)^2 + (σ3 − σ1)^2)/ 2) [/math]

Donde σ1, σ2 y σ3 son los autovalores de σ (también conocidos como tensiones principales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico. Se utiliza en ingenería para predecir el punto de fallo de un material cuando está sujeto a cargas complejas. La fórmula de la tensión de Von Mises ayuda a sumar todas las diferentes formas de fuerza que afectan al material y ofrece un número que dice cuánto puede soportar el material antes de que se rompa. 

h = 2/10;  % Paso de muestreo
ro = 1:h:2;  % Asignación de las variables rho y theta
zeta = linspace(atan(1/2), pi-atan(1/2), pi/h);
[RR, Z] = meshgrid(ro, zeta);  % Mallado
lro=length(ro);
lzeta=length(zeta);
%creamos las componentes de la matriz sigma
a=@(RR,Z) (1./RR).*exp(RR-1).*cos(2.*Z-(pi/2)); %Elemento (1,1) y (3,3) de la matriz
b=@(RR,Z) (1/2).*exp(RR-1).*sin(2.*Z-(pi/2)); %Elemento (1,2) y (2,1) de la matriz
c=@(RR,Z) 3.*(1./RR).*exp(RR-1).*cos(2.*Z-(pi/2)); %Elemento (2,2) de la matriz
Msig=[];
%hallamos los autovalores
for i=1:lzeta
 for j=1:lro
    A=a(RR(i,j),Z(i,j));
    B=b(RR(i,j),Z(i,j));
    C=c(RR(i,j),Z(i,j));
    Msig=[A,B,0;B,C,0;0,0,A];
    autovalores=eig(Msig);
   VM=sqrt(((autovalores(1)-autovalores(2))^2+((autovalores(2)- ...
       autovalores(3))^2+((autovalores(3)-autovalores(1))^2))*1/2));
   G(i,j)=VM;
 end
end
x=RR.*cos(Z);  % Cambio a coordenadas cilíndricas
y=RR.*sin(Z);
%dibujamos la tensión de Von Mises
subplot(1,2,1)
surf(x,y,G);
axis equal
xlabel('x')
ylabel('y')
view(2)
axis([-3,3,-1,3])
title('Tensión de Von Mises')
colorbar
shading interp
subplot(1,2,2)
surf(x,y,G);
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
title('Tensión de Von Mises')
colorbar
shading interp
%valor máximo de la tensión de Von Mises
max=max(max(G));
disp(max)


izquierda centro

El valor máximo de la tensión de Von Mises es de 3.3291

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Análisis_y_visualización_de_transformaciones_en_un_cuarto_de_anillo&oldid=86362»