Campos y temperatura en anillo plano (Grupo 14A)

Trabajo realizado por estudiantes
Título	Campos y temperatura en anillo plano (Grupo 14A)
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Jose Manuel Alonso de Caso, Raquel Álvarez,Alejandro del Rey, Raphael Gide,Miguel López, Alejandro Mendizabal.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 . Mallado del anillo plano

Partimos de una superficie plana circular centrada en el origen, delimitada por los radios de valores 1 y 2 unidades, que expresaremos en coordenadas cilíndricas para el estudio en este trabajo.

Mallado del anillo plano. Paso de muestreo h=0.1

%vectores de las coordenadas
rho=linspace(1,5,50);
phi=linspace(0,2*pi,50);


%malla
h=0.1;     % Paso de muestreo
%coordenadas polares
r=1:h:2;                            
t=0:h:2*pi+h;          
[rr,tt]=meshgrid(r,t);   %mallado      
%Parametrización
x=rr.*cos(tt);                 
y=rr.*sin(tt);
mesh(x,y,0*x)           
axis([-3,3,-3,3])

2 . Distribución de temperaturas en el anillo

La temperatura se distribuye dentro del anillo siguiendo la función [math]T=(8-y^2+2y)exp(-x^2+2x-1)[/math]. La temperatura en cada punto del disco debido a esta distribución se puede intuir mediante las curvas de nivel de la representación gráfica adjunta abajo. Se trata de una distribución cuyo foco de calor se encuentra desplazado con respecto al origen de coordenadas. De esta forma se observan curvas concéntricas que están formadas por puntos que se encuentran a la misma temperatura (curvas de nivel).

Utilizaremos la expresión [math] max(max(T)) [/math] para obtener el valor máximo de la temperatura (T) para que matlab nos lo muestre. El máximo de temperatura en el dominio es [math] 8.9945 [/math], aproximadamente [math]9[/math]

h=0.1;                               
r=1:h:2;                            
t=0:h:2*pi+h;          
[rr,tt]=meshgrid(r,t); %mallado    
x=rr.*cos(tt);                 
y=rr.*sin(tt);
mesh(x,y,0*x);
axis([-3,3,-3,3])    
T=(8-y.^2+2.*y).*exp(-x.^2+2.*x-1);   % Expresión del campo escalar térmico
subplot(2,1,1)
surf(x,y,T)          
axis([-3,3,-3,3])  
subplot(2,1,2)
contour(x,y,T,25);
axis([-3,3,-3,3])

Distribución de temperaturas y curvas de nivel del anillo.

3 . Gradiente de la temperatura

El gradiente de la temperatura [math]∇T=(T_1,T_2)= \frac{\partial T}{\partial x} \vec i + \frac{\partial T}{\partial y} \vec j = (-2x+2).exp(-x^2+2x-1) \vec i + (-2*y+2).exp(-x^2+2*x-1) \vec j [/math] se define como ortongonal a las lineas de nivel, como podemos comprobar visualmente en la figura. Se observa también que el gradiente en cada punto converge hacia el foco de calor.

Gradiente de la temperatura en el anillo

h=0.1;
r=1:h:2;
t=0:h:2*pi+h;
[r,t]=meshgrid(r,t);
x=r.*cos(tt);
y=r.*sin(tt);
%Relación polares-cartesianas
u=sqrt(x.^2+y.^2);
v=atan(y./x); 
f=(8-y.^2+2.*y).*exp(-x.^2+2.*x-1);
fx=(-2.*x+2).*exp(-x.^2+2.*x-1);     % Derivada parcial en x
fy=(-2.*y+2).*exp(-x.^2+2.*x-1);     % Derivada parcial en y
hold on
quiver(x,y,fx,fy)
axis([-3,3,-3,3]) 
contour(x,y,f,25)    %Curvas de nivel
axis([-3,3,-3,3])
hold off

4 . Desplazamiento provocado por [math]\vec u[/math]

A continuación aplicamos el desplazamiento provocado por el vector [math]\vec u = (1-ρ)^2 g_ρ [/math] en el mallado y lo representamos gráficamente. Cada vector representa la fuerza aplicada en cada punto del sólido y nos indica la dirección de desplazamiento que lleva cada punto.

Desplazamiento provocado por el vector u en el anillo.

h=0.1;
r=1:h:2;
t=0:h:2*pi+h;
[rr,tt]=meshgrid(r,t);
x=rr.*cos(tt);
y=rr.*sin(tt);
u=sqrt(x.^2+y.^2);
v=atan(y./x); 
ux=((1-u).^2).*x./u;  
uy=((1-u).^2).*y./u; 
quiver(x,y,ux,uy)
axis([-3,3,-3,3])

5 . Sólido antes y después del desplazamiento

Obtenemos ahora la posición final de los puntos del anillo cuando actúa en él el desplazamiento anterior. Comparando el estado inicial y el estado final del anillo, vemos que las partículas del sólido se mueven en la dirección de los vectores del campo u, como habíamos visto en el apartado anterior. El desplazamiento provoca que el anillo se expanda hasta llegar al radio exterior de valor 3, siendo anteriormente 2. El radio interior se mantiene en 1.

%malla
h=0.1;    %paso de muestreo 
%coordenadas polares      
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t); %mallado
x=rr.*cos(tt);       
y=rr.*sin(tt);
u=sqrt(x.^2+y.^2);
v=atan(y./x); 
%Desplazamientos originados por la vibración
ux=((1-u).^2).*x./u;  
uy=((1-u).^2).*y./u; 
%Posición final de los puntos de la placa
dx=x+ux;
dy=y+uy;
subplot(1,2,1) %antes
mesh(x,y,0.*x)       
axis([-3,3,-3,3])      
subplot(1,2,2)  %despues
mesh(dx,dy,0.*dx)
axis([-3,3,-3,3])

expansión de la malla al aplicar [math] \vec u [/math]

6 . Divergencia del campo vectorial [math]\vec u[/math]

Para analizar el cambio de volumen debido al desplazamiento anterior utilizamos el operador diferencial de la divergencia \( \nabla \cdot \vec u \): [math] ∇·\vec u=\frac{ 1}{ √g}·\frac{\partial }{\partial x^i}·(√g·u^i) = [/math]

Representación de la divergencia del campo u

Representación de la divergencia de u en 3D

h=0.1;                               
r=1:h:2;                            
t=0:h:2*pi+h;          
[rr,tt]=meshgrid(r,t);     
x=rr.*cos(tt);                 
y=rr.*sin(tt);
mesh(x,y,0*x)
axis([-3,3,-3,3])    
div=1./rr-4+3.*rr;   % Expresión de la divergencia
surf(x,y,div)          
axis([-3,3,-3,3])

7 . Rotacional [math]\vec u[/math]

Para aplicar la formula del rotacional necesitamos las covariantes del campo [math] \vec u⃗ = (1-ρ)^2 . g^ρ = u_i . g^i[/math] que al aplicar la matriz de gram de las coordenadas polares queda [math] \vec u⃗ = (1-ρ)^2 . g_ρ = u^i . g_i[/math]

[math] \nabla\times \vec u⃗ =\frac{ 1}{ ρ}\left| \begin{matrix} \ g_ρ & \ g_θ & \ g_z \\ & & \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ & & \\ u_1 & u_2 & u_3 \end{matrix}\right| [/math]

[math] ∇\times\vec u=(\frac{\partial u^3}{\partial ρ}-\frac{\partial u^2}{\partial θ})·g_ρ-(\frac{\partial u^3}{\partial ρ}-\frac{\partial u^1}{\partial z})·g_θ+(\frac{\partial u^2}{\partial ρ} -\frac{\partial u^1}{\partial θ})·g_z [/math] Las coordenadas en [math]g_ρ[/math] y [math]g_θ[/math] s anulan al derivar, y en la [math]g_z[/math]:

[math] [ -(1-\frac{ 1}{ ρ^2})*senθ + (1-\frac{ 1}{ ρ^2})*senθ ]·g_z [/math] [math] ∇φ=0·g_i [/math]

8 . Obtencion del campo de tensiones [math]σ[/math]

Calculamos el tensor gradiente del campo [math]\vec u[/math] y a continuación su parte simétrica que llamaremos [math]ε[/math]. Definimos el tensor de tensiones [math]σ[/math] como [math]σ = λ∇·u1+2µε[/math] siendo los coeficientes de Lamé [math] λ [/math] y [math] µ[/math] iguales a [math]1[/math].

Representación de las tensiones normales en dirección gρ.

h=0.1;                         
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t);     
x=rr.*cos(tt);     %Relación polares-cartesianas  
y=rr.*sin(tt);
sigma1= 1./rr-8+7.*rr; % Campo de tensiones en dirección gρ
surf(x,y,sigma1)              
axis([-3,3,-3,3])          
axis equal
colorbar