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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad (Grupo 9A)
Asignatura	Teoría de Campos
Curso	2016-17
Autores	Yasmine El Khattabi, Luis Rincón Bagüés, Javier Balairón Garcia, David Carracedo Esteban
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 .Introducción

El siguiente trabajo versa sobre la aplicación en un cuarto de anillo circular de una placa un campo de temperaturas y otro de desplazamientos. La placa ocupa un cuarto de anillo circular comprendido entre los radios 1 y 3 del primer cuadrante. El campo U es un campo en coordenadas cilidricas de la siguiente forma

u = f(ρ)∙g_θ

y el campo u calculado es:

u (ρ,θ) = 0.1(ρ-1) 

La función de temperaturas es f=log(y+2). Además se han estudiado la divergencia del campo, el rotacional, las tensiones de deformaciones, la tensión de Von Mises y se ha realizado el calculo de la masa de la placa.

2 .Mallado

En primer lugar se define la placa (anillo comprendido entre las circunferencias de radios 1 y 3 en el primer cuadrante). Para ello se utilizan coordenadas cilíndricas con rho perteneciente al intervalo [1,3] y theta perteneciente a [0,pi/2]. Es importante que el muestreo sea múltiplo de pi/2 para que la placa quede definida completamente.

h=(pi/2)/20;
u=1:h:3;
v=0:h:pi/2;
[rho,theta]=meshgrid(u,v);
xx=rho.*cos(theta);
yy=rho.*sin(theta);
mesh(xx,yy,0*xx) 
axis([-1,4,-1,4])

3 .Temperatura

La temperatura sigue la ecuación f=log(y+2). Para representarla utilizamos un mallado con el comando meshgrid.

h=0.1;
x=-1:h:4;
y=x;
[mx,my]=meshgrid(x,y);
f=log(my+2);
mesh(mx,my,f);

4 .Gradiente

Para representar el gradiente utilizamos el comando gradient. Para comprobar que es perpendicular al campo de temperaturas lo representamos en planta con el comando contour y colocamos juntos ambos gráficos con el comando subplot.

h=0.1;
x=-1:h:4;
y=x;
[mx,my]=meshgrid(x,y);
f=log(my+2);
[px,py]=gradient(f);
hold on
contour(mx,my,f);
quiver(mx,my,px,py);
hold off

5 .Campo de desplazamientos U.

Para la obtención del campo de desplazamientos u, se han tenido en cuenta las condiciones propuestas por el enunciado:

Los puntos situados en ρ = 1 no sufren desplazamiento.
|∇xu| = (3ρ-2)/10 
u = f(ρ)∙g_θ

Resolviendo analiticamente se ha llegado al siguiente resultado:

u (ρ,θ) = 0.1(ρ-1) 

6 .Campo de vectores en los puntos del mallado del sólido.

Partiendo del campo de desplazamientos, vamos a representar el campo de vectores con este codigo de matlab Nota: gθ = -ρ·sen(θ)i + ρ·cos(θ)j

clear all 
h=0.1;
u = 1:h:3;
v = 0:h:pi/2;
[rho,theta] = meshgrid(u,v);
xx = rho.*cos(theta);
yy = rho.*sin(theta);
U = 0.1.*(rho-1);
U1 = U.*(-yy);
U2 = U.*(xx);
mesh(xx,yy,0.*xx)
view(2)
axis([-1,4,-1,4])  
hold on
quiver(xx,yy,U1,U2)
title('campo de vectores u')

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

h=0.1;                  
u=1:h:3
v=0:h:pi/2+h;          
[rho,theta]=meshgrid(u,v); 
xx=rho.*cos(theta);        
yy=rho.*sin(theta);
U = 0.1.*(rho-1);
U1 = U.*(-yy);
U2 = U.*(xx);
Dx=xx+U1;
Dy=yy+U2;
subplot(1,2,1);
mesh(xx,yy,0.*xx)
view(2);
axis([-1,4,-1,4]);
title('Antes de desplazamiento')
subplot(1,2,2);
mesh(Dx,Dy,0.*Dx)
view(2);   
title('Despues de desplazamiento')
axis([-1,4,-1,4]);

8 .Representación de la divergencia de U

El campo es el siguiente: u (ρ,θ) = 0.1(ρ-1)

Al calcular la divergencia, el resultado es 0 debido a que el campo u=f(ρ)∙g_θ es una función de ρ y al derivar en el calculo de la divergencia es nulo.Por tanto,la placa no sufre incrementos de volumen sino que se deplaza permaneciendo su volumen constante.

9 .Representación del rotacional de U.

Partiendo del campo: u (ρ,θ) = 0.1(ρ-1) y calculando analiticamente obtenemos el rotacional: (3ρ-2)/10. En representación con Matlab se utiliza el siguiente código:

rot=(3*rho-2)/10;
h=(pi/2)/20;
u=1:h:3;
v=0:h:pi/2;
[rho,theta]=meshgrid(u,v);
xx=rho.*cos(theta);
yy=rho.*sin(theta);
mesh(xx,yy,rot)

Se observa que los puntos que sufren un mayor rotacional son los con ρ=3.

10 . Tensor de deformaciones.

11 .Tensiones ortogonales.

12 .Tensión de Von mises.

La tensión de Von Mises viene dada por la siguiente expresión: [\sigma_{VM} = \sqrt{\frac{(\sigma_1-\sigma_2)^2+ (\sigma_2-\sigma_3)^2+ (\sigma_3-\sigma_1)^2}{2}}\]. Se tomarán como autovalores de la matriz de tensiones los siguientes, correspondientes a los valores de sigma_1, sigma_2 y sigma_3: sigma_1=0; sigma_2=1/10; sigma_3=-1/10. Estos autovalores se han obtenido partir de la matriz del tensor de tensiones. Sustituyedo en la fórmula de Von Mises obtenemos el valor de sqrt(3)/10. Por lo tanto la tensión es un escalar y afecta a la placa por igual en todos sus puntos.

13 . Masa de la placa.

Para el cálculo de la masa total del anillo, se emplea la función densidad [math]d(x,y)=1+e(-abs(x)/(y+1)^2)[/math] que se integra a lo largo del dominio para obtener de esa manera una aproximación a dicha magnitud. El código de Matlab empleado es el siguiente:

N1=200; 
N2=100;
a=1; 
b=3; 
c=0; 
d=pi/2;
h1=(b-a)/N1; h2=(d-c)/N2;
u=a:h1:b; 
v=c:h2:d;
[rho,theta]=meshgrid(u,v);
f=1+exp(-abs(rho)./(theta+1).^2);
w1=ones(N1+1,1); 
w(1)=1/2; w(N1+1,1)=1/2;
w2=ones(N2+1,1);
w(1)=1/2; w(N2+1,1)=1/2;
result=h1*h2*w2'*f*w1

El resultado es 4.8128 unidades de masa.