Visualización de campos escalares y vectoriales en elasticidad Grupo 21C

1 Introducción

Este estudio consistirá en analizar el comportamiento de una placa circular frente a la acción de los diferentes campos que actúan sobre ella. Para ello los analizaremos y estudiaremos observando su variación mediante la representación gráfica de estos con la ayuda del programa informático MATLAB.

2 Placa de estudio

Consideramos una placa plana (en dimensión 2) que ocupa el anillo comprendido entre las circunferencias centradas en el origen de radios 2 y 1. Dibujamos un mallado que representa los puntos interiores del sólido tomando como paso de muestreo h=1/10 para las variables x e y.

Placa Circular

h=0.1;
 u=1:h:2;
 v=0:h:2*pi+h;
 %mallado de la gráfica
 [uu,vv]=meshgrid(u,v);
 figure(1)
 xx=uu.*cos(vv);
 yy=uu.*sin(vv);
 mesh(xx,yy,0*xx)
 axis([-3,3,-3,3])
 view(2)

3 Campo de Temperatura

En ella tenemos definida la temperatura, que proviene de un foco de calor muy concentrado en el origen y viene dada en función de ρ y θ como T(ρ,θ)=−log(ρ+0.1). Observamos en la gráfica su concentración en el origen, lo que hace que tome sus valores máximos en ρ=0. La función esta expresada en polares y para trabajar en Matlab hay que pasarlo a coordenadas cartesianas. T(x,y)=-log⁡(√(x^2+y^2 )+0.1)

Gráfica de la Temperatura

%Relación polares-cartesianas
 r=sqrt(xx.^2+yy.^2);
 tt=atan(yy./xx);
 %Campo escalar T
 figure(2)
 T=-log(0.1+r);
 surf(xx,yy,T)
 axis([-3,3,-3,3])
 view(2)

3.1 Variación de temperatura

Calculamos el gradiente de T para observar su variación y representamos también sus curvas de nivel, observando que estas son perpendiculares a los vectores que representan el campo vectorial de ∇T.

La función T está expresada en polares. Las curvas de nivel, junto con el gradiente, quedarían:

Gradiente y Curvas de Nivel

Tx=-1./(r.*(0.1+r)).*xx;
 % derivada parcial en x
 Ty=-1./(r.*(0.1+r)).*yy;
 % derivada parcial en y
 subplot(1,2,1),quiver(xx,yy,Tx,Ty)
 axis([-3,3,-3,3])
 % Región del gráfico
 subplot(1,2,2),contour(xx,yy,T,10)
 %dibujo curvas de nivel
 axis([-3,3,-3,3])
 view(2)
 % representación campo vectorial u
 figure(4)
 % componente y de u
 ux=sin(pi.*tt./2)./(30.*r).*(xx./r);
 % componente y de u
 uy=sin(pi.*tt./2)./(30.*r).*(yy./r);
 quiver(xx,yy,ux,uy)
 axis([-3,3,-3,3])

4 Campo de vibraciones

Consideramos ahora el campo de vectores [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}. [/math] que determina el desplazamiento que sufre cada punto del sólido cuando este actúa sobre la placa.

Dibujamos a continuación el campo de vectores en los puntos del mallado del sólido.

Campo Vectorial

figure(9)
 m=inline('((sin((pi*atan(yy./xx)/2)))./(30*(xx.^2+yy.^2))).*xx','xx','yy');
 n=inline('((sin((pi*atan(yy./xx)/2)))./(30*(xx.^2+yy.^2))).*yy','xx','yy');
 M=m(xx,yy);
 N=n(xx,yy);
 quiver(xx,yy,M,N);
 view(9);

Una vez este campo de fuerzas actúa sobre el sólido, éste sufre un desplazamiento tal y como observamos a continuación:

Desplazamiento de la malla

u=1:0.1:2;
 v=0:0.1:2*pi+0.1;
 [uu,vv]=meshgrid(u,v);
 figure(1) 
 subplot(1,2,1);
 xx=uu.*cos(vv);
 yy=uu.*sin(vv);
 mesh(xx,yy,0*xx)
 axis([-3,3,-3,3])
 view(2) 
 subplot(1,2,2);
 ux=((sin(pi*vv/2))./(30*uu)).*cos(vv);
 uy=((sin(pi*vv/2))./(30*uu)).*sin(vv);
 mesh(xx+ux,yy+uy,0*xx);
 axis([-3,3,-3,3]);

4.1 Divergencia de [math]\vec u[/math]

Para calcular la divergencia vamos a usar las coordenadas polares a través de la expresión: [math]\nabla \cdot \vec u = \frac{1}{√g}\frac{∂}{∂x^i}(√gu^i)[/math]

La divergencia es nula. Por tanto podemos decir que el campo no es ni una fuente ni un sumidero y además el cuerpo no sufre ningún cambio de volumen debido al desplazamiento.

4.2 Rotacional de [math]\vec u[/math]

Al calcular el rotacional del campo [math]|\nabla \times \vec u|[/math] observamos que los puntos con mayor rotacional son aquellos en los que la malla ha sufrido el desplazamiento.

Rotacional del Campo u

figure(12);
 %Función que resulta al multiplicar vectorialmente el operador Nabla y el campo u
 r=inline('(pi*(cos(pi/2*atan(y./x))))./(60*sqrt(x.^2+y.^2))','x','y');
 R=r(xx,yy);
 surf(xx,yy,R);
 axis([-3,3,-3,3]);
 view(12);

5 Tensiones

El campo vectorial de tensiones se reparte por la placa tal y como observamos en la siguiente figura, donde se muestra que ésta está concentrada en tres focos principales donde las tensiones son mayores.

Campo de tensiones

5.1 Módulo de Tensiones

Los valores que adoptan las tensiones en cada punto se pueden observar representando su módulo en una gráfica. De esta forma vemos claramente los tres focos principales en 'rojo' mientras que en el borde de la placa, donde esta toma unos colores 'azulados' las tensiones son casi nulas.

Módulo de las tensiones

Módulo de tensiones visto en planta

h=0.1;
 u=1:h:2;
 v=0:h:(2*pi+h);
 [uu,vv]=meshgrid(u,v);
 xx=uu.*cos(vv);
 yy=uu.*sin(vv);
 ww= -(1./(15*(uu.^2))).*(sin((pi*vv)/2)).*(cos(vv))-(pi./(60*(uu))).*(cos((pi*vv)/2)).*(sin(vv));
 zz=-(1./(15*(uu.^2))).*(sin((pi*vv)/2)).*(sin(vv))+(pi./(60*(uu))).*(cos((pi*vv)/2)).*(cos(vv));
 figure(1)
 quiver(xx,yy,ww,zz)
 axis([-3,3,-3,3])
 M=sqrt(((1./(15*(uu.^2))).*(sin((pi*vv)/2))).^2+(((pi./(60*(uu))).*(cos((pi*vv)/2))).^2)./(uu.^2));
 figure(2)
 surf(xx,yy,M)
 figure(3)
 surf(xx,yy,M)
 axis([-3,3,-3,3])