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 circular plana centrada en el origen, delimitada por las circunferencias con radios de 1 y 2 unidades. Utilizaremos un mallado que expresaremos en coordenadas polares, sean ρ є [1,2] y θ є [0,2Π) para el estudio en este trabajo.

Mallado del anillo plano. Paso de muestreo h=0.1

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)e^{-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 sendas curvas concéntricas que están formadas por puntos que se encuentran a la misma temperatura (curvas de nivel).

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

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

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]) 
 % Expresión del campo escalar térmico   
T=(8-y.^2+2.*y).*exp(-x.^2+2.*x-1);   
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])

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).e^{-x^2+2x-1} \vec i + (-2*y+2).e^{-x^2+2*x-1} \vec j [/math] se define como ortogonal 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. La representación del campo vectorial gradiente indica, en cada punto del sólido, la dirección en la cual la temperatura aumenta más rápido, y su módulo nos dirá la rápidez con la que aumenta la temperatura en esa dirección.

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 uno.

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 [math] \vec u [/math], 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.

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

%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])

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] La divergencia es mayor en el borde exterior del anillo del dominio. En los puntos del limite interior del anillo no hay desplazamiento con respecto a la posición inicial y podemos ver que ésta se incrementa progresivamente en los puntos que están más hacia el exterior. Nuestra figura no tiene volumen puesto que es una superficie plana.

derecha

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].

El valor de la matriz del tensor es :[math]\begin{bmatrix} \frac{ 1}{ ρ}}-8+7ρ & 0\\ 0 & \frac{ 3}{ ρ}}-8+5ρ\\ \end{bmatrix} [/math]

sabiendo que el tensor se puede expresar como: [math]σ= (frac{ 1}{ ρ}}-8+7ρ).g_ρ\otimes g_ρ + (\frac{ 3}{ ρ}}-8+5ρ).g_θ\otimes g_θ [/math]

8.1 . Tensiones normales en dirección g_ρ

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

8.2 . Tensiones normales en dirección g_θ/ρ

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);     
figure(1)
x=rr.*cos(tt);     %Relación polares-cartesianas  
y=rr.*sin(tt);
sigma2=3./rr-8+5.*rr; % Campo de tensiones en dirección gθ/ρ
surf(x,y,sigma2)              
axis([-3,3,-3,3])         
axis equal
colorbar

Conclusión: la tensión es más elevada en la direccion g_ρ que en la g_θ por tanto es mayor en la dirección radial que en las perpendiculares a los radios.

9 . Tensiones tangenciales

Las tensiones tangenciales son cero porque la matriz del tensor de tensiones,σ, es diagonal. Los elementos asociados a la dirección tangencial de la tensión, es decir, los que no pertenecen a la diagonal principal, son cero, no existen en las direcciones de g_θ y g_ρ. Por tanto solo podemos observar las tensiones normales.

[math]\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho = 0 [/math]

10 . Tensión de Von Mises

[math]\sigma_{VM} = \sqrt{\frac{(\sigma_1-\sigma_2)^2+ (\sigma_2-\sigma_3)^2+ (\sigma_3-\sigma_1)^2}{2}}[/math]

Siendo [math]\sigma_1=\frac{ 1}{ ρ}}-8+7ρ, \sigma_2=\frac{ 3}{ ρ}}-8+5ρ, y \sigma_3=0,[/math] que por adaptación al plano escribiremos como: [math]\sigma_{VM} = \sqrt{\frac{(\sigma_1-\sigma_2)^2}{2}}[/math]

que al sustituir los autovalores, queda:

[math]\sigma_{VM} = \frac{(-2+2.ρ^2)}{\sqrt{2}.ρ}[/math]

Representación de las tensiones de Von Mises

h=0.1;                         
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t);     
figure(1)
x=rr.*cos(tt);     
y=rr.*sin(tt);
fvm=(-2+2.*rr.^2)./(rr.*sqrt(2));
surf(x,y,fvm)              
axis([-3,3,-3,3])         
axis equal
colorbar

11 . Masa del anillo

Utilizamos el siguiente código de Matlab para aproximar la integral:

Representación de la función de masa en cada punto.

N1=200; N2=100;                  %Numero de puntos
a=1; b=2; c=1; d=2;             %Extremos del intervalo
h1=(b-a)/N1; h2=(d-c)/N2;
u=a:h1:b; v=c:h2:d;              %coordenadas iniciales
[uu,vv]=meshgrid(u,v);           %coordenadas del rectángulo
f=vv.*uu.*exp(-1./uu.^2);                %función
w1=ones(N1+1,1);                 %vector de masas
w(1)=1/2; w(N1+1)=1/2;
w2=ones(N2+1,1);                 %vector de masas
w(1)=1/2; w(N2+1)=1/2;
result=h1*h2*w2'*f*w1

La masa del anillo con función de densidad [math] d(x,y,z) = xye−1/x2 \lt\math\gt es 1.4635[/math]