Visualización de Campos en un Anillo Plano (Grupo15-C)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo 15-C
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Jaime Peña, Iñigo Uraga, Pablo Molinero, Iñigo Diez, Daniel Pacheco y Pablo Vazquez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 . Introducción

Mallada anillo plano

Dada una placa plana con forma de anillo circular centrado en el origen, y comprendido entre los radios 1 y 2 (Figura 1), se estudiarán las variaciones de temperatura en este objeto y las deformaciones causadas por el campo:

[math] \vec u(\rho,\theta)={(1-\rho)^2}\vec g_{\rho} [/math]

En esta placa circular consideramos dos magnitudes físicas: la temperatura [math]T(\rho,\theta,t)[/math], que depende de las variables espaciales [math]ρ[/math] y [math]θ[/math] y del tiempo [math]t[/math], y los desplazamientos [math] \vec u (ρ,θ,t)[/math], que dependen de los mismos parámetros que la temperatura. De esta forma, si definimos [math] \vec r_0(\rho,\theta)[/math] como el vector de posición de los puntos de la placa en reposo, la posición de cada punto [math](\rho,\theta)[/math] de la placa para un instante de tiempo [math]t[/math] viene dada por:

[math] \vec r (\rho,\theta,t)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta,t) [/math]

2 . Mallado del anillo plano

Tomamos una superficie circular plana centrada en el origen, comprendido por las circunferencias de radios 1 y 2 unidades. Utilizaremos un mallado que expresaremos en coordenadas polares, sean ρ є [1,2] y θ є [0,2Π).

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

3 . Distribución de temperaturas en el anillo

La temperatura del sólido viene dada por la función: [math]T=(8-y^2+2y)e^{-x^2+2x-1}[/math] La representación de las curvas de nivel de la temperatura está representada en la gráfica inferior. El foco de calor viene representado en rojo y se observa que no esta en el origén. De esta forma se observa que la temperatura se propaga con curvas concéntricas isotermas.

Utilizando Matlab obtenemos que el máximo de la temperatura 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]) 
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])

4 . 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] es ortogonal a las lineas de nivel. El gradiente indica, en cada punto del sólido, la dirección en la cual la temperatura aumenta más rápido, es decir hacia el foco, 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);
u=sqrt(x.^2+y.^2); %Relación polares-cartesianas
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

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

A continuación aplicamos el campo de desplazamientos provocado por el vector :[math] \vec u(\rho,\theta)={(1-\rho)^2}\vec g_{\rho} [/math] en el mallado y lo representamos gráficamente.

Desplazamiento en el anillo provocado por el vector u

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

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

Representamos la posición final de los puntos del anillo cuando actúa en él el desplazamiento anterior. Observamos que las partículas del sólido se mueven en la dirección de los vectores del campo [math] \vec u [/math]. El desplazamiento hace que el anillo se agrande llegando a un radio exterior de 3 unidades, siendo anteriormente 2 unidades, manteniendo el radio interior en 1 unidades.

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

h=0.1;    %Paso de muestreo 
r=1:h:2;   %Coordenadas polares             
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);
ux=((1-u).^2).*x./u;  %Desplazamientos originados por la vibración
uy=((1-u).^2).*y./u; 
dx=x+ux; %Posición final de los puntos de la placa
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])

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

Calculamos la divergencia de u en todos los puntos del sólido, utilizando el operador \( \nabla \cdot \vec u \): [math] ∇·\vec u=\frac{ 1}{ √g}·\frac{\partial }{\partial x^i}·(√g·u^i) [/math] Los puntos de mayor divergencia están en la frontera exterior del anillo. En los puntos de la frontera interior no hay desplazamiento y observamos que aumenta progresivamente conforme se acercan a la frontera exterior. No se aprecia cambio de volumen, ya que se trata de una representación de una superficie plana.

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

8 . 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} \\ & & \\ (1-ρ)^2 & 0 & 0 \end{matrix}\right| =0 [/math]

9 . Obtención del campo de tensiones [math]σ[/math]

Calculamos el tensor gradiente del campo [math]\vec u[/math] mediante la derivada covariante (forma para obtener el gradiente de un campo vectorial)y a continuación obtenemos 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] Ambas componentes representan [math]σ[/math] en dirección [math]g_ρ[/math] y [math]g_θ[/math], respectivamente.

9.1 . Tensiones normales en dirección [math]g_ρ[/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

9.2 . Tensiones normales en dirección [math]g_θ/ρ[/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);     
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.

10 . Tensiones tangenciales

Las tensiones tangenciales son cero ya que la matriz de tensiones,[math] σ [/math], es diagonal debido a que el campo de desplazamiento produce una expansión en la dirección estrictamente radial sin producir ninguna torsión alrededor del eje vertical del origen de coordenadas. 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 tensiones en las direcciones de los planos ortogonales a g_θ y g_ρ.

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

11 . Tensión de Von Mises

La tensión de Von Mises es determinada por:

[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ρ, \sigma_3=0,[/math] los autovalores de la matriz de tensiones. Adaptando la fórmula a nuestro caso, anillo circular plano, la desarrollamos como: [math]\sigma_{VM} = \sqrt{\frac{(\sigma_1-\sigma_2)^2}{2}}[/math]

Sustituyendo los autovalores:

[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

12 . Masa del anillo

Obtenemos la masa del anillo utilizamos la función de densidad [math] d(x,y,z) = xye^{-1/x^2} [/math] que integramos en todo dominio para hallar una aproximación de esta magnitud.

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 en el primer cuadrante es 1.4635. Debido a su simetría impar, si integraramos en todo el dominio la masa total sería 0. Como dos de los cuadrantes son de masa idéntica pero negativa, al multiplicar por 4 la masa del primer cuadrante, obtenemos lo que sería una suma de la masa en valor absoluto de todo el dominio. [math] 1.4635 \times 4 = 5.854 = Masa Total [/math]