Visualización de campos escalares y vectoriales en elasticidad. (Grupo 45)

Trabajo realizado por estudiantes
Título	Campos de temperaturas y deformaciones en un cuarto de anillo (Grupo 45)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Jose Andres Bello Amado Javier Diaz Santamaria Pelayo Gomez Lobo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 . Mallado de la placa

Tenemos que generar el mallado del cuarto de anillo circular en Matlab, para esto emplearemos el comando meshgrig. Anillo está centrado en (0,0) y sus radios son R1 y R2 que valen 1 y 2 respectivamente. Nos pide que la representación este encuadrada en la región \((x, y) ∈ [-3; 3] × [−1; 3]\).

El paso de muestreo debe ser de \(h =\frac{2}{10}\)

Pasamos a coordenadas cilíndricas estas coordenadas cartesianas, ρ es el radio del anillo \(θ\) y el ángulo que forman las rectas límite del cuarto de anillo con los ejes, de valor:

Las rectas límite son: \(y =\frac{x}{2}\) y \(y =\frac{-x}{2}\)

\(α = tg(\frac{1}{2}) ≈ 26,6°\)

Por lo tanto las coordenadas en cilíndricas son:

[math](ρ,θ) ∈ [1; 2] × [α,π-α][/math]

Siendo α=26,6°.

Figura 1

clear
clc

% Variables
h = 2/10;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

x = D.*cos(A);
y = D.*sin(A);
mesh(x,y,0.*x)
view(2)
axis([-3,3,-1,3])

title('Cuarto de Anillo Circular')

2 . Curvas de temperatura

La distribución de temperatura en este sólido queda definida por el campo escalar T:

[math]T(x,y)=cos(x^2)+sin((y−1)^2)[/math]

Ahora podemos definir el gradiente de T como T. En este caso, indica dirección en la que aumenta temperatura, y el modulo del gradiente indica cuanto lo hace. Formula del gradiente en cartesianas:

[math]\nabla T = \frac{\partial T}{\partial x}\vec i+\frac{\partial T}{\partial y}\vec j+\frac{\partial T}{\partial z}\vec k[/math]

Con lo que llegamos a:

[math]\nabla T = -2xsin(x^2)\vec i+2cos((y-1)^2)(y-1)\vec j[/math]

clear
clc

% Placa
h = 0.2;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

X = D.*cos(A);
Y = D.*sin(A);

% Campo de temparatura
T = cos(X.^2)+sin((Y-1).^2);

% Curvas de nivel
figure
mesh(X,Y,0*X);
hold on
contour(X,Y,T,20)
colorbar
title('Curvas de temperatura')
axis ([-3, 3, -1, 3,])
view(0,90);
hold off
Tmax = max(max(T));

Se observa que las temperaturas más altas están en la parte superior del sólido, siendo esta de 1.8415 Ahora podemos comprobar como el gradiente es ortogonal a las curvas de nivel:

clear
clc

% Placa
h = 0.2;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

x = D.*cos(A);
y = D.*sin(A);

% Campo de temparatura
T = cos(x.^2)+sin((y-1).^2);

% Curvas de nivel
figure
mesh(x,y,0*x);
hold on

% Gradiente
X = -2.*x.*sin(x.^2);

3 . Ley de Fourier

La conducción térmica sigue la Ley de Fourier, dice que el flujo de trasferencia de calor es el contrario al gradiente en esa dirección. Imaginemos una placa sólida con superficie A situada entre dos grandes placas separadas por una distancia Y. Supongamos que antes de t=0, el material tiene una temperatura To. En t=0, la placa experimenta un calentamiento repentino a la temperatura Ts, que se mantiene constante. Con el tiempo, el perfil de temperatura cambia en relación con la posición hasta alcanzar un perfil lineal de temperaturas (estado estacionario) independiente del tiempo y la posición después de un período de transición. Después de este cambio, se requiere un flujo de calor constante para mantener las temperaturas. Esta condición se cumple para gradientes de temperatura bajos, aplicaciones de manera diferencial en tres dimensiones (coordenadas rectangulares).

[math]\vec Q =-k▽T \left\{\begin{matrix} \vec Q = conductividad térmica del material \\\\\\\\ k = constate de conductividad térmica del material \end{matrix}\right.[/math]

En este caso suponemos que k=1, por lo tanto, la fórmula de la energía calorífica queda:

[math]\vec Q=-▽T[/math]

Lo único que quedaría por hacer es cambiar el signo del gradiente como muestra la ecuación. Lo hacemos a continuación y obtenemos el campo vectorial de la energía calorífica:

[math]\vec Q =2xsen(x^2)\vec i -2(y-1)cos((y-1)^2)\vec j [/math]

clear
clc

% Variables
h = 2/10;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

% Conversión a cartesianas
x = D.*cos(A);
y = D.*sin(A);
mesh(x,y,0*x);

% Aplicación 
X = 2.*x.*sin(x.^2);
Y = 2.*(y-1).*cos((y-1).^2);
hold on
quiver(x,y,X,Y);
hold off
view(2)
axis([-3,3,-1,3])

title('Energía Calorífica')

4 . Campo de vectores en mallado solido

Debemos representar el campo de vectores deslizantes en nuestro cuarto de anillo, para ello nos dan un campo vectorial [math] \vec u = \frac{1}{2}sen(θ) \vec e_ρ [/math] .

Pasamos [math]\vec u[/math] a cartesianas, sabiendo que [math] \vec e_ρ= cos(θ) \vec i + sin(θ)\vec j[/math]

clear
clc

% Variables 
h = 2/10;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

x = D.*cos(A);
y = D.*sin(A);
mesh(x,y,0.*x)

view(2)
axis([-3,3,-1,3])
hold on
X=0.5.*sin(A).*cos(A);
Y=0.5.*sin(A).*sin(A);
quiver(x,y,X,Y);

title('Campo de Vectores')

5 . Imagen solido antes y después desplazamiento

Pasamos a representar el sólido sin ser deformado y también una vez deformado, dicha deformación la produce el campo de presiones [math] \vec u [/math] en cartesianas. Este campo [math] \vec u [/math] es el que hemos calculado en el apartado 4. Como [math] \vec e_{\theta} [/math] es nula, la deformación se producirá en la dirección de [math] \vec e_{\rho} [/math].

clear
clc

% Sin Deformar
h = 2/10;                         
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);
x = D.*cos(A);
y = D.*sin(A);
subplot(1,2,1)                  
V = mesh(x,y,0*x);
view(2)
set(V,'EdgeColor','y');
axis ([-3,3,-1,3])
title('Placa Sin Deformar');

% Deformada
subplot(1,2,2)                     
N = (sin(A).*cos(A))./2;
M = (sin(A).^2)./2;
X = x+N;
Y = y+M;
K = mesh(X,Y,0*X);
view(2)
set(K,'EdgeColor','r');
axis ([-3,3,-1,3])
title('Placa Deformada');

% Comparación
hold on
figure                                 
K = mesh(X,Y,0*X);
axis ([-3,3,-1,3])
view(2)
set(K,'EdgeColor','r');
title('Comparación de Ambas Placas');
hold on
V = mesh(x,y,0*x);
set(V,'EdgeColor','y');
hold off

6 . Divergencia

Calculamos la divergencia a partir de la siguiente fórmula:

[math]\nabla ·u = \frac{1}{\rho} [\frac{\partial }{\partial \rho} (\rho u_\rho)+\frac{\partial }{\partial θ} (u_θ)+\frac{\partial }{\partial z} (\rho u_z)][/math]

y de esto, operando, llegamos a:

[math]\nabla ·u = \frac{senθ}{2\rho}[/math]

A partir de esto, representamos en Matlab la divergencia gráficamente:

clear
clc

% Variables
h = 0.2;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A]=meshgrid(dist,alpha);

X=D.*cos(A);
Y=D.*sin(A);

% Divergencia
divergencia=(sin(A)./2./D);
surf(X,Y,divergencia)
view(0,90);
axis equal
axis([-3,3,-1,3])
colorbar
title('Divergencia')

% Máximos y Mínimos
Divergenciamax=max(max(divergencia));
Divergenciamin=min(min(divergencia));

6.1 Máximos y mínimos

A partir del gráfico, observamos que las divergencias mayores están en la parte central del sólido y la mínima en los dos extremos:

.Divergencia máxima: 0,5000

.Divergencia mínima: 0,1118

7 . Rotacional del campo

clear
clc

% Variables
h = 2/10;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha); 

X = D.*cos(A);
Y = D.*sin(A);
mesh(X,Y,0*X);
 
% Rotacional
rotacional = 1/2.*cos(A)./D;
 
% Grafica
surf(X,Y,rotacional);
view(2)
axis([-3,3,-1,3])
colorbar

title('Rotacional del Campo')

Observamos, que los puntos de mayor rotacional son aquellos con menor ángulo [math]\theta[/math] y menor valor de ρ .

8 . Tensiones normales

clear
clc

%Variables 
h = 2/10;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

X=D.*cos(A);
Y=D.*sin(A); 

Z=(1./(2.*D)).*(sin(A));
subplot(3,1,1);
surf(X,Y,Z);
axis equal
title('Tensión normal en ρ'); 

W=((1/2).*sin(A))+sin(A);
subplot(3,1,2);
surf(X,Y,W);
axis equal;
title('Tensión normal en θ'); 

M=(1./(2.*D)).*(sin(A));
subplot(3,1,3);
surf(X,Y,M);
axis equal
title('Tensión normal en z');

9 . Tensiones tangenciales

clear
clc

%Variables
h = 2/10;
dist = 1:h:2;
w = atan(1/2);
alpha = linspace(w,pi-w,pi/h);
[D,A] = meshgrid(dist,alpha);

X=D.*cos(A);
Y=D.*sin(A); 

% Tensiones 
M = cos(A)./(2.*D);
surf(X,Y,M);
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
title('Tensiones Tangenciales');

10 . Tensión de Von Mises

right

clear
clc
h=2/10;  
dist=1:h:2;  
w = atan(1/2);
alpha=linspace(w,pi-w,pi/h);
[D,A]=meshgrid(dist,alpha);  
x=D.*cos(A);  
y=D.*sin(A);

sig= [];  
V=zeros(length(alpha),length(dist));
a = @(R,T) (sin(R)./(2.*T));
b = @(R,T) (cos(R)./(2.*T));
c = @(R,T) ((3.*sin(R))./(2.*T));

for i=1:length(alpha)  
   for j=1:length(dist)
     sig(1,1)= a(D(i,j),A(i,j));
     sig(1,2)= b(D(i,j),A(i,j));
     sig(1,3)= 0;
     sig(2,1)= b(D(i,j),A(i,j));
     sig(2,2)= c(D(i,j),A(i,j));
     sig(2,3)= 0;
     sig(3,1)= 0;
     sig(3,2)= 0;
     sig(3,3)= a(D(i,j),A(i,j));
     
     [v,u]=eig(sig);  
     V(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);  
   end
end
subplot(1,3,1)  
surf(x,y,V)
shading interp
view(2)
axis equal
colorbar;
title('XOY')

subplot(1,3,2)  
surf(x,y,V)
shading interp
axis equal
axis vis3d
colorbar;
title('Tensión de Von Mises')

subplot(1,3,3)  
hold on
m=V(1:size(V,1),1);
theta=x(1:length(x),3);
surf(x,y,V)
view(0,0)
shading interp
MAX = max(m);

for k = 1:length(m)  
   if m(k) == MAX
     plot3(theta(k),0,MAX,'xr','markersize',10)  
   end 
end 
axis equal
colorbar
title('XOZ')
hold off