Visualización de campos vectoriales y escalares en elasticidad. Grupo C-16

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad.Grupo C-16.
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Iván Díez Berjano (1111) Alejandro Montesinos Pérez (954) Jorge Martinez Velasco (931) Pablo Sobrado (976) Luis Vergara Rojo (1031) Ignacio Perez Abellan (928)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa circular plana que ocupa el anillo circular centrado en el origen y comprendido entre los radios 1 y 2. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura (T(x,y)=e^-y) y los desplazamientos (ū(ρ,θ)=(1-ρ)²·g_θ), dependientes cada una de las variables que se mostraran más adelante en los apartados correspondientes.

1 - Mallado de los puntos interiores del sólido

Mallado de la corona circular.

El mallado consiste en la representación de los puntos que componen la placa sólida en forma de corona circular. Tomaremos como paso de muestreo 0.1 para las variables X e Y.

%Definición de parámetros de la superficie
u=[1:0.1:2];          
v=[0:0.1:2*pi+0.1];
%Creación del mallado
[uu,vv]=meshgrid(u,v); 
%Parametrización
xx=uu.*cos(vv);        
yy=uu.*sin(vv);
%Dibujo de la malla
mesh(xx,yy,0*xx)  
%Ajustamos la gráfica 
axis([-4,4,-4,4]); 
view(2);

La parametrización quedará hecha únicamente en este apartado, no se escribirán en el resto de programas.

2 - Campo de temperaturas del sólido

La temperatura del sólido varía a lo largo de este según la función T(x,y)=e^-y. Dibujaremos las curvas de nivel para observar en que lugares esta temperatura es máxima. La función de temperaturas depende solo de la variable Y, por lo que se mantendrá constante para todo valor de X.

Distribución de la temperatura en la placa

%Campo escalar temperatura
T=exp(-yy); 
%Creación de la primera subventana
subplot(1,2,1);
%Representación de la malla del campo de temperaturas
mesh(xx,yy,T);
%Vista según los planos X e Y    
view(2)
axis([-4,4,-4,4]) 
%Creación de la segunda subventana
subplot(1,2,2);
%Representación de la misma malla
%Vista en 3D
mesh(xx,yy,T); 
colorbar;

En la gráfica observamos que cuanto mayor es el valor de la Y, menor es la temperatura; luego los puntos de máxima temperatura son los más alejados del foco y por tanto con menor valor de Y.
LA TEMPERATURA SERÁ MÁXIMA EN EL PUNTO (X,Y,Z)=(0,-2,7.5)

3 - Gradiente del campo de temperaturas

El gradiente del campo de temperaturas (∇T) nos indica cómo varía la temperatura al trasladarnos por el plano. El módulo del gradiente nos indica la rapidez con la que aumenta la temperatura. El gradiente apunta en la dirección en la cual la temperatura aumenta más rápido, que es la dirección en donde la derivada direccional será máxima. Como hemos visto en teoría, el vector gradiente será perpendicular a las líneas de contorno de nuestro sólido (Podemos verlo en la figura mostrada a la derecha). En la representación podemos observar las diferentes lineas circulares y concéntricas que muestran la disminución de la temperatura a medida que los puntos del sólido se alejan del foco.

Representación del campo vectorial del gradiente de T (∇T)

%Gradiente del campo
Gradiente=-exp(-yy);	
%Creación de la primera subventana
subplot(1,2,1)
%Representación del gradiente junto con el campo de temperaturashold on
%Malla del campo de temperaturas
mesh(xx,yy,T);
%Gradiente
quiver3(xx,yy,T,0*xx,Gradiente,0*T);
hold off
%Creación de la segunda subventana
subplot(1,2,2)
%Representación del gradiente 
quiver(xx,yy,0*xx,Gradiente); colorbar

4 - Campo de desplazamientos ū

Usaremos como campo de desplazamientos ū(ρ,θ)=(1-ρ)²·g_θ. Dibujaremos su correspondiente campo de vectores. Sabemos que g_θ = -ρ·sen(θ)i + ρ·cos(θ)j. Por ello divididermos y definimos el campo en función de los dos componente i,j. Definiremos pues, fx y fy respectivamente.

Desplazamiento del campo vectorial ū(ρ,θ)=(1-ρ)²·g_θ

%Componente x(uu,vv) de mi campo vectorial
fx=(1-uu).^2.*(-(uu.*sin(vv))); 
%Componente y(uu,vv) del campo vectorial
fy=(1-uu).^2.*(uu.*cos(vv));  
%Representamos el campo vectorial sobre la placa
quiver(xx,yy,fx,fy)  
%Límites de los ejes x e y en la representación de la placa
axis([-3,3,-3,3])
%Visualización
view(2)         
%Estos apartados siguientes podríamos omitirlos, pero como los empleamos a la hora de hacer el programa y la correspondiente captura 
%de pantalla en matlab, los pondremos a modo explicativo, pero no son necesarios.
%Nombramos el eje x
xlabel('x');   
%Nombramos el eje y
ylabel('y');  
%Damos un título al gráfico
title('Mallado y campo de vectores u')

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

Representaremos en una misma ventana nuestro sólido antes y después del desplazamiento dado por el campo de vectores anterior y así poder compararlos. Definiremos las funciones de desplazamiento (desx, desy) en función de ambas componentes, sumando las componentes del vector de desplazamientos a las coordenadas respectivas de posición. (como vemos en las líneas 8 y 10, de nuestro programa)

Placa del sólido antes y después desplazamiento

%Creación de la primera subventana
subplot(1,2,1)
%Representación
mesh(xx,yy,0*xx)      
axis([-4,4,-4,4]);
view(2);
%Desplazamiento en eje x
desx=xx+fx;
%Desplazamiento en eje y
desy=yy+fy;
%Creación de la segunda subventana
subplot(1,2,2)
%Representación del desplazamiento
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); view(2);

Podemos apreciar que el desplazamiento es mínimo.

6 - Divergencia

Aquí calculamos la divergencia del campo de vectores u en todos los puntos del sólido y la representamos. A través de esta gráfica veremos qué puntos tienen mayor divergencia.

Divergencia del campo del desplazamineto

%Calculamos la divergencia a mano y el resultado fue cero
%Dibujamos la superficie
div=0.*yy
surf(xx,yy,div)
%Visualización en 2D
view(2)
colorbar

La divergencia del campo de vectores u es cero, por lo que en la gráfica no se aprecia variación. Considerando la divergencia como una medida del cambio de volumen local debido al desplazamiento, al ser esta nula, no habrá ningunn cambio de volumen a causa del desplazamiento.

7 - Rotacional

En este apartado calculamos el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. A través de la gráfica se observa cuáles son los puntos que sufren un mayor rotacional. Hemos calculado el rotacional a mano, en un papel aparte, obteniendo que es Rot [ ū(ρ,θ) ] = 4ρ² - 6ρ + 2 .

Rotacional del campo del desplazamineto

%Creamos la matriz del rotacional 
rotacional = zeros(64,11);
%Asignamos los valores del rotacional a la matriz componente a componente
%Calculamos así el módulo del rotacional
i=1;
j=1;
for i = 1:64
  for j = 1:11
    rotacional(i,j)= 2+4.*(uu(i,j).^2)-6.*uu(i,j);
    j= j+1;
  end
  i = i+1;
end
%Creamos la primera subventana
subplot(1,2,1)
%Representamos el campo vectorial
quiver3(xx,yy,xx*0,xx*0,xx*0,-rotacional); 
%Creamos la segunda subventana
subplot(1,2,2)
%Dibujamos el mallado del rotacional
mesh(xx,yy,rotacional); colorbar;

Los puntos que sufren un mayor rotacional son aquellos que mayor coordenada Z tienen.

8 - Tensor de tensiones

Definimos la parte simétrica del tensor gradiente de ū y calculamos el tensor de tensiones σ. Dibujamos las tensiones normales en la dirección que marca los vectores i y j, los cuales en nuestro trabajo sustiruiremos por g_ρ y g_θ/ρ respectivamente, y comparamos esta gráfica con las del módulo de la divergencia y el rotacional.

Tensiones normales, Divergencia y Rotacional

subplot(1,3,1)
%Modulo de las tensiones normales
tensN=0.*uu;
%Dibujamos la malla
surf(xx,yy,tensN)                 
subplot(1,3,2)
%Modulo de la divergencia
div=0.*uu;
%Dibujamos la superficie del divergente                     
surf(xx,yy,div);          
subplot(1,3,3)
%Malla del rotacional
mesh(xx,yy,rotacional)
colorbar

9 - Tensiones tangenciales respecto al plano ortogonal a g_ρ

Aquí dibujamos las tensiones tangenciales respecto al plano ortogonal a g_ρ, junto a la divergencia y el rotacional. A modo de explicación podríamos decir que las tensiones tangenciales representadas hacen referencia a esfuerzos cortantes del material.

Comparacion de las tensiones tangenciales con las del módulo de la divergencia y las del módulo del rotacional

subplot(1,3,1)
%Modulo de las tensiones tangenciales
tensT=uu.*(1-uu).*(-2+uu.*(1-uu).^2);                    
%Dibuja la malla de las tensiones tangenciales
mesh(xx,yy,tensT) 
%Vista desde el plano ortogonal
view(2)
subplot(1,3,2)
%Modulo de la divergencia
div=0.*uu;  
%Malla de la diverencia
mesh(xx,yy,div);          
subplot(1,3,3)
%Malla del rotacional
mesh(xx,yy,rotacional)
colorbar

En la gráfica apreciamos las zonas de mayor tensión tangencial. Esta zona es casi la más externa de la corona circular, sin llegar a ser lo más externa posible pues como vemos en el dibujo empezaría a disminuir, tomando un color amarillo, en vez de rojo.

10 - Tensiones tangenciales respecto al plano ortogonal a g_θ/ρ

Representamos las tensiones tangenciales por el plano ortogonal a g_θ/ρ, al igual que en el apartado anterior.

Comparacion de las tensiones tangenciales con las del módulo de la divergencia y las del módulo del rotacional

subplot(1,3,1)
%Modulo de las tensiones tangenciales
tensT=uu.*(1-uu).*(-2+uu.*(1-uu).^2);                    
%Dibujamos la malla de las tensiones tangenciales
surf(xx,yy,tensT)
%Hemos variado la posición de la figura para verla como se nos pide.
subplot(1,3,2)
%Tenemos el módulo de la divergencia
div=0.*uu;  
%Dibujamos la malla de la divergencia
mesh(xx,yy,div);          
subplot(1,3,3)
%Dibujamos la malla del rotacional
mesh(xx,yy,rotacional)
colorbar

En la gráfica apreciamos donde las tensiones son mayores. A mayor altura (z) estas tensiones serán también mayores.

11 - Tensión de Von Mises

Calculamos la tensión de Von Mises a partir de la fórmula dada calculando los autovalores de σ y la representamos. Los autovalores calculados son los de la matriz de tensiones. Los hemos calculado aparte en matlab usando la sugerencia del comando 'Eig'.

Representación de la Tensión de Von Mises.

%Calculamos la tension de Von Mises
VM=sqrt(3).*(1-uu).*(-2+uu.*(1-uu).^2);
%Dibujamos la superficie
surf(xx,yy,VM)
colorbar

En la gráfica podemos apreciar en que punto la tensión es mayor.

12 - Masa total

A partir de la función de densidad d(x,y,z)=xylog(x+2) calculamos la masa y la representamos en el gráfico de la derecha, que aún no siendo necesario nos proporciona una idea un poco más visual. Con el método del trapecio hallamos una aproximación del valor de la integral, dividiendo la proyección del sólido en cuadrados de lado h=1/1000 para mayor precisión.

Distribución de la masa de la placa

%En primer lugar vamos a dibujar la gráfica de distribución de la masa
%Masa=densidad*Volumen
%Sabemos que densidad=xx*yy*log(xx+2);
%Calculamos el volumen a mano
masa=xx.*yy.*log(xx+2).*(pi*4-pi*1);               
%Representación de la malla
surf(xx,yy,masa)                     
axis([-2,2,-2,2])                   
colorbar 

%En segundo lugar, vamos a calcular la integral de la densidad
%Método del trapecio
%Definimos las particiones de los rectángulos
h=1/1000;
%Definimos nuestra función de densidad
f=(uu.*vv.*exp((-1)./(uu.^2)));
a=(h^2)*f;
%Nos muestra el resultado
masa=sum(sum(a))

LA MASA TOTAL ES: M=0.0021