Estudio de las transformaciones de un cuarto de anillo (Grupo 6)

Trabajo realizado por estudiantes
Título	Estudio de las transformaciones de un cuarto de anillo (Grupo 6)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Ismael García Sanz Daniel Pérez Cámara Mohammed Rami Alejandro Brea Gilardi
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El trabajo 6 consiste en realizar lo mismo que en el trabajo 1 pero tomando 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. Dibujar con los ejes en el cuadrado [−3, 3] × [−1, 3] .

Tomaremos el campo [math] \vec u(ρ,θ) = \frac{1}{2}sin(θ)\vec e_p [/math]

Tomaremos una función de temperatura [math]T(x,y)=cos(x^2)+sen((y-1)^2)[/math]

En este caso y para simplificar la mayoría de los cálculos se trabaja en coordenadas cilíndricas.

A continuación tenemos un índice donde se puede ver el contenido del trabajo

1 . Mallado de la placa

El anillo circular viene definido por dos parámetros: estar comprendido entre los radios 1 y 2 y y en en plano y ≥|x|/2. Debido a que es un anillo circular, será mas sencillo trabajar en coordenadas cilíndricas.

Del primer parámetro obtenemos el valor de ρ, que se encuentra entre [1,2]. Del segundo parámetro obtenemos el valor de θ, que se encuentra entre [arctg(1/2), π - arctg(1/2)], que numéricamente sería [0.463647609,pi-0.463647609].

Por lo que en coordeneadas cilíndrica se representa como (ρ,θ)∈ [1,2] × [0.463647609,π-0.463647609].

Representación del anillo circular

clear
close all 
%Apartado 1
%Parametrizacion con Muestreo
h=2/10;
%Dimensiones del anillo circular
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Mallado del anillo circular
[R,T]=meshgrid(r,t);
%Cambio de Coordenadas Cartesianas a Coordenadas Cilindricas
X=R.*cos(T);
Y=R.*sin(T);
Z=Y*0;
%Representación del anillo circular
mesh(X,Y,Z);
%Longitud y nombre de los ejes
axis equal
axis([-3,3,-1,3]); 
title("Cuarto de Anillo Circular") 
xlabel('Eje X')
ylabel('Eje Y')
%Vista en 2D
view(2);

2 . Curvas de temperatura

La temperatura del sólido viene dada por el siguiente campo escalar:

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

2.1 Gráfico representación de la temperatura y de las curvas de nivel en T

En las siguientes gráficas se ven representada la función temperatura en el anillo circular y las curvas de nivel en T.

Las zonas azuladas son las zonas donde la temperatura es menor, mientras que en las zonas amarillentas la temperatura es mayor. Las zonas blancas es donde la temperatura es nula.

Representación de las curvas de nivel

clear
close all 
%Apartado 2
%Parametrizacion con Muestreo
h=2/10;
%Dimensiones del anillo circular
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Mallado del anillo circular
[R,T]=meshgrid(r,t);
%Cambio de Coordenadas Cartesianas a Coordenadas Cilindricas
X=R.*cos(T);
Y=R.*sin(T);
Z=Y*0;
%Funcion de Temperatura
F=cos(X.^2)+sin((Y-1).^2);
%Representación el campo escalar de temperaturas
hold on
subplot(1,2,1)
surf(X,Y,F)  
view(2)
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
colorbar
%Gráficas de curvas de nivel
subplot(1,2,2)
contour(X,Y,F,20/h)
title('Curvas de Nivel')
colorbar
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
hold off

2.2 Temperatura Máxima

La temperatura máxima es 1,8415ºC, que corresponde al valor 83 de F, siendo X(83)=1,22246e - 16 y Y(83)=2.

max(max(T))
ans=1,8415

2.3 Gradiente de T

El gradiente de la temperatura representa la dirección de más velocidad de cambio de la función.

El gradiente se define como:

[math]\nabla T(x,y) =\frac{∂T}{∂X}\vec i + \frac{∂T}{∂Y}\vec j [/math]

Por lo tanto, el gradiente será:

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

Se sabe que el gradiente siempre es ortogonal a las curvas de nivel, como se ve tanto en 2D como en 3D.

Gradiente de T

clear
close all 
%Gradiente
h=2/10;  %Paso de Muestreo
r=1:h:2;  %Asignación de las variables rho y theta
t=linspace(0.463647609,pi-0.463647609,pi/h);
[R,T]= meshgrid(r,t);  %Creación del mallado
X=R.*cos(T);  %Cambio a coordenadas cilíndricas
Y=R.*sin(T);
Z=Y*0;
%Gradiente de la temperatura en x
dx=-2.*X.*sin(X.^2); 
%Gradiente de la temperatura en y
dy=2*cos((Y-1).^2).*(Y-1); 
%Función de temperatura
F=cos(X.^2)+sin((Y-1).^2); 
%Gráfica en 3D
subplot(1,2,2)  
hold on
surf(X,Y,F)
quiver3(X,Y,F,dx,dy,0*F)
axis vis3d
view(3)
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Gradiente 3D')
hold off
%Gráfica en 2D
subplot(1,2,1)  
view(2)
hold on
contour(X,Y,F,60)
quiver(X,Y,dx,dy)
axis([-3,3,-1,3])
title('Gradiente de la temperatura en 2D')
hold off

3 . Ley de Fourier

Le energía calorífica es la energía que posee un cuerpo en forma de calor.

De acuerdo a la Ley de Fourier la energía calorífica [math]\vec Q~[/math] viaja de acuerdo a la fórmula:

[math]\vec Q=-k∇T[/math]

donde k es la constante de conductividad térmica de la placa que supondremos κ = 1.

La energía calorífica [math]\vec Q~[/math] en este caso es:

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

Debido a que la energía calorífica depende de un gradiente, conserva la propiedad de que es ortogonal a las curvas de nivel.

Fourier

clear
close all 
%Apartado 3
%Parametrizacion con muestreo h=2/10
k=1
h=2/10;
%Dimensiones del anillo circular
r=1:h:2;  
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Mallado del anillo circular
[R,T]= meshgrid(r,t);
%Cambio de Coordenadas Cartesianas a Coordenadas Cilindricas  
X=R.*cos(T);  
Y=R.*sin(T);
Z=Y*0;
%Función de temperatura
F=cos(X.^2)+sin((Y-1).^2); 
contour(X,Y,F,20);
%Gradiente de T
%Gradiente de la temperatura en x
dx=-2.*X.*sin(X.^2); 
%Gradiente de la temperatura en y
dy=2*cos((Y-1).^2).*(Y-1); 
%Calculo de energia calorifica
Qx=-k*dx;
Qy=-k*dy;
%Titulo y dimension de ejes
title('Campo Vectorial Q');
xlabel('Eje X');
ylabel('Eje Y');
%Barra de color dependiendo de la temperatura
hold on
quiver(X,Y,Qx,Qy)
axis equal
colorbar

4 . Campo de vectores en mallado solido

Dibujar el campo de vectores en los puntos del mallado del sólido, en [math]t=0[/math].

Se hará una representación del campo de vectores mediante el vector: [math]\vec u(ρ,θ) = \frac{1}{2}sin(θ)\vec e_ρ[/math].

Definiremos los intervalos de nuestras variables y después generaremos el mallado para obtener las matrices [math]ro[/math] y [math]theta[/math]. Estableceremos las variables [math]X[/math] e [math]Y[/math] y mediante ellas definiremos el campo de desplazamientos con las variables [math]FX[/math] y [math]FY[/math]. Después representaremos el campo de vectores de desplazamiento ajustando los ejes según los pedidos.

Representación del campo de vectores

%Se generan las variables
h=0.2;
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Se hace el mallado 
[ro,theta]=meshgrid(r,t);
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(1/2).*sin(theta).*cos(theta);
FY=(1/2).*sin(theta).*sin(theta);
%visualizacion del dibujo
quiver(X,Y,FX,FY);
axis equal
axis([-3,3,-1,3]);
view(2)
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos sobre la curva')

5 . Imagen solido antes y después desplazamiento

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec{u}[/math] (en [math]t=0[/math]). Dibujar ambos en la misma figura usando el comando [math]subplot[/math].

Usaremos el vector: [math]\vec u(ρ,θ) = \frac{1}{2}sin(θ)\vec e_ρ[/math] y empezaremos definiendo la región en la que se encuentra el sólido y después calcularemos el sólido antes y después del desplazamiento. Para terminar realizaremos una comparación del movimiento del sólido.

%Se generan las variables
h=0.2;
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Se hace el mallado 
[ro,theta]=meshgrid(r,t);
X=ro.*cos(theta);
Y=ro.*sin(theta);
FX=(1/2).*sin(theta).*cos(theta);
FY=(1/2).*sin(theta).*sin(theta);
%Sólido antes del desplazamiento
figure
surf(X,Y,0*X);
view(2)
axis equal
axis([-3,3,-1,3]);
title('Antes de desplazamiento');
xlabel('Eje x');
ylabel('Eje y');
%Sólido después del desplazamiento
figure
UX=X+FX;
UY=Y+FY;
surf(UX,UY,0*UX);
view(2)
axis equal
axis([-3,3,-1,3]);
title('Después de desplazamiento');
xlabel('Eje x');
ylabel('Eje y');
%Comparación
figure
hold on
a=mesh(X,Y,0*X);
b=mesh(UX,UY,0*UX);
set(a,'Edgecolor','r');
set(b,'Edgecolor','b');
view(2)
axis equal
axis([-3,3,-1,3]);
title('Comparación');
xlabel('Eje x');
ylabel('Eje y');
legend('Sin desplazamiento','Con desplazamiento');
hold off

Sólido antes del desplazamiento

Sólido después del desplazamiento

Comparación

6 . Divergencia

Dibujar [math]∇·\vec{u}[/math] en [math]t=0[/math]. Determinar analíticamente los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. La divergencia es una medida del cambio de volumen local debido al desplazamiento. ¿Se puede apreciar esto en la gráfica?

Calculamos la divergencia mediante la fórmula: [math]∇·\vec{u} = \frac{1}{ρ}[\frac{∂}{∂ρ}(p \vec u_ρ)+\frac{∂}{∂θ}(\vec u_θ)+\frac{∂}{∂z}(p \vec u_z)][/math].

Comenzaremos definiendo los intervalos de las variables, generaremos el mallado y calcularemos la divergencia mediante la formula de antes. Después dibujaremos la divergencia y calcularemos la divergencia máxima y mínima siendo el valor máximo 0.5 y el valor mínimo 0.1118.

Divergencia del campo vectorial

%Se generan las variables
h=0.2;
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Se hace el mallado 
[ro,theta]=meshgrid(r,t);
X=ro.*cos(theta);
Y=ro.*sin(theta);
%Divergencia
div=(1./ro).*(1/2.*sin(theta));
%Dibujo
surf(X,Y,div);
view(2)
axis equal
axis([-3,3,-1,3]);
xlabel('Eje x');
ylabel('Eje y');
colorbar
%Divergencias máximas y mínimas
DivMax=max(div);
DivMin=min(div);
title('Divergencia campo vectorial')

7 . Rotacional del campo

Calcular [math]|∇ × \vec{u}|[/math] en todos los puntos del sólido en [math]t = 0[/math] y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

El rotacional se haría con la siguiente fórmula: [math]∇ × \vec{u}=\frac{1}{ρ}\begin{vmatrix}\vec{e_ρ} & ρ\vec{e_θ}& \vec{e_z}\\ \frac{\partial }{\partial ρ} & \frac{\partial }{\partial θ}& \frac{\partial }{\partial z}\\ u_{ρ}& ρu_{θ}& u_{z} \end{vmatrix} = \frac{1}{ρ}\begin{vmatrix} \vec{e_ρ} & ρ\vec{e_θ} & \vec{e_z}\\ \frac{\partial }{\partial ρ} & \frac{\partial }{\partial θ} &\frac{\partial }{\partial z} \\ 1/2sinθ & 0 & 0\end{vmatrix}[/math]

[math] \mathbb{\nabla \times \vec u} = {1 \over ρ} \; \begin{vmatrix} {\vec e_ρ} & {ρ· \vec e_θ} & {\vec e_z}\\ {\partial \over \partial ρ} & {\partial \over \partial θ} & {\partial \over \partial z}\\ {\frac{1}{2}sen(θ)} & {0} & {0}\\ \end{vmatrix} [/math] = [math] -{1 \over ρ}{\partial \over \partial θ}{\frac{1}{2}sen(θ)}{\vec e_z} = \frac{-1}{2ρ} \cos(θ){\vec e_z}[/math]

El resultado de esa operación es: [math] |∇ ×\vec u |[/math]= [math]\frac{1}{2ρ}\cos(θ) [/math]

Rotacional del campo de vectores

clear
close all 
%Parametrizacion con muestreo h=2/10
k=1;
h=2/10;
%Dimensiones del anillo circular
r=1:h:2;  
t=linspace(0.463647609,pi-0.463647609,pi/h);
%Se crea el mallado
[ro,theta]=meshgrid(r,t);
x=ro.*cos(theta);
y=ro.*sin(theta);
%El rotacional ya calculado
rotacional=(1/2).*(cos(theta)./ro);
%Gráfica en 2D
subplot(1,2,1);
surf(x,y,rotacional);
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);
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)

El valor del rotacional máximo alcanzado es de: 0.4472

8 . Tensiones normales

Definamos [math]ϵ(\vec{u}) = (∇\vec{y} + ∇\vec{u}^t)/2[/math], la parte simétrica del tensor gradiente de [math]\vec{u}[/math] conocido como tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]σij[/math] a través de la fórmula [math]σ = λ∇ · \vec{u} 1 + 2µϵ[/math] donde [math]1[/math] es el tensor identidad en el conjunto de vectores libres del espacio [math]R^3[/math] y [math]λ[/math], [math]µ[/math] son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. A pesar de que los desplazamientos son planos (es decir [math]\vec{u}[/math] no tiene componente en la dirección de [math]\vec{k}[/math]) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa. Tomando [math]λ = µ = 1[/math], dibujar las tensiones normales en la dirección que marca el eje [math]\vec{i}[/math], es decir [math]\vec{i}· σ ·\vec{i}[/math], las tensiones normales en la dirección que marca el eje [math]\vec{j}[/math], es decir [math]\vec{j} · σ · \vec{j}[/math] y las correspondientes al eje [math]\vec{k}[/math], es decir [math]\vec{k}· σ · \vec{k}[/math] (dibujar las que no son nulas).

Se realizarán las tensiones normales según las diferentes direcciones [math]\vec e_ρ[/math], [math]\vec e_θ[/math] y [math]\vec e_z[/math].

[math]ε=\begin{pmatrix} \ 0 & \frac{cos(θ)}{4ρ} & 0 \\ \frac{cos(θ)}{4ρ} & \frac{sen(θ)}{2ρ} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

La divergencia es [math]\nabla\cdot\vec u=\frac{sen(θ)}{2ρ}[/math]

Entonces σ es igual a : [math] σ= \begin{pmatrix} \frac{sen(θ)}{2ρ} & 0 & 0 \\ 0 & \frac{sen(θ)}{2ρ} & 0 \\ 0 & 0 & \frac{sen(θ)}{2ρ} \end{pmatrix} + 2\begin{pmatrix} \ 0 & \frac{cos(θ)}{4ρ} & 0 \\ \frac{cos(θ)}{4ρ} & \frac{sen(θ)}{2ρ} & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} \ \frac{sen(θ)}{2ρ} & \frac{cos(θ)}{2ρ} & 0 \\ \frac{cos(θ)}{2ρ} & \frac{3sen(θ)}{2ρ} & 0 \\ 0 & 0 & \frac{sen(θ)}{2ρ} \end{pmatrix}[/math]

Para [math]\vec e_ρ[/math] habrá que hacer: [math]\vec e_ρ[/math] × σ × [math]\vec e_ρ[/math] = [math]\frac{sen(θ)}{2ρ}[/math]

Para [math]\vec e_θ[/math] habrá que hacer: [math]\vec e_θ[/math] × σ × [math]\vec e_θ[/math] = [math]\frac{3sen(θ)}{2ρ}[/math]

Para [math]\vec e_z[/math] habrá que hacer: [math]\vec e_z[/math] × σ × [math]\vec e_z[/math] = [math]\frac{sen(θ)}{2ρ}[/math]

%la matriz del gradiente es el siguiente 
%A=[0,cos(theta)/2,0;0,0,0;0,0,0];
%el tensor de deformación es 
%E=(A+A')/2;
%tensor de tensiones
%T=[sin(theta)/2*ro,0,0;0,sin(theta)/2*ro,0;0,0,sin(theta)/2*ro]+2.*E;
%Tensiones normales en la dirección del eje e? 
r=1:0.1:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
[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=T.*[1;0;0];%elemento(1,1,1);
a=sin(theta)./(2.*ro);
hold on 
figure(1);
surf(x,y,a);
axis equal
title('Tensión normal en la dirección e_?')
view(3);
colorbar;
%Tensiones normales en la dirección del eje e? ?
%se calculan las componentes de la matriz sigma
b=(3/2).*(sin(theta)./ro);
figure(2);
surf(x,y,b);
axis equal
title('Tensión normal en la dirección e_?')
view(3);
colorbar;
%Tensiones normales en la dirección del eje e? z
%se calculan las componentes de la matriz sigma
c=sin(theta)./(2.*ro);;%elemento(3,3,3)
figure(3)
surf(x,y,c);
axis equal
title('Tensión normal en la dirección e_k')
view(3);
colorbar;
hold off

9 . Tensiones tangenciales

Calcular las tensiones tangenciales respecto al plano ortogonal a [math]\vec e_ρ[/math]. Se pide además las no nulas.

[math]|σ·\vec e_ρ-(\vec e_ρ·σ·\vec e_ρ)·\vec e_ρ| = \left | \begin{pmatrix} \ \frac{sen(θ)}{2ρ} & \frac{cos(θ)}{2ρ} & 0 \\ \frac{cos(θ)}{2ρ} & \frac{3sen(θ)}{2ρ} & 0 \\ 0 & 0 & \frac{sen(θ)}{2ρ} \end{pmatrix}·\begin{pmatrix} 1\\0\\0 \end{pmatrix} - \frac{sen(θ)}{2ρ}·\begin{pmatrix} 1\\0\\0 \end{pmatrix}\right | = \left | \begin{pmatrix} 0\\ \frac{cos(θ)}{2ρ}\\0 \end{pmatrix}\right| [/math]

Tensiones tangenciales

% Se realiza el mallado
h=0.1;
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
[rr,tt]=meshgrid(r,t); 
xx=rr.*cos(tt);
yy=rr.*sin(tt);
% Tensión tangencial
tension=cos(tt)/2;
% Tensión tangencial en 2D
subplot(1,2,1)
surf(xx,yy,tension)
axis equal
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
view(3)
colorbar 
title('Tensión tangencial en 3D en la dirección del eje e_?');

10 . Tensión de Von Mises

La tension de Von Mises es una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro). Está definida por la siguiente fórmula:

[math]\ σ_{VM}= \sqrt {\frac{(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}{2}}[/math],

donde [math]\ σ_1,σ_2,σ_3, [/math] son los autovalores del tensor de tensiones [math]σ(ρ,θ)[/math]

La tensión máxima que puede soportar el cuarto de anillo antes de deformarse es de 0.7746.

h=0.2;
r=1:h:2;
t=linspace(0.463647609,pi-0.463647609,pi/h);
[RR,TT]=meshgrid(r,t);
xx=RR.*cos(TT);
yy=RR.*sin(TT);
zz=zeros(size(xx));

M11 = @(R,T) sin(T)./(2.*R);
M12 = @(R,T) cos(T)./2;
M22 = @(R,T) sin(T).*R./2;

%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) = M12(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
figure;
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');

figure
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

figure

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

derecha

izquierda

centro