Visualización de campos escalares y vectoriales en elasticidad (Grupo C4)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad (Grupo C4)
Asignatura	Teoría de Campos
Curso	2017-18
Autores	Pablo Salvadores, Raquel Marco
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Para la visualización de campos escalares y vectoriales, consideraremos una placa definida por el anillo circular centrado en el origen y comprendidos entre los radios 1/2 y 2, en el plano y≥0. Con esta definición, tendremos una placa formada por un semianillo circular. Consideraremos ademas como ejes los pertenecientes al rectángulo [-3,3]x[-1,3]. Para simplificar la tarea de definir la placa, usaremos coordenadas cilindricas.

En la placa vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x,y)[/math] que depende de las dos variables espaciales (x,y), y los desplazamientos [math]u(x,y)[/math] producidos por la acción de una fuerza determinada. De esta forma, si definimos r0(x,y) como el vector posición de los puntos de la placa antes de la deformación, la posición de cada punto (x,y) de la placa después de la deformación viene dada por : [math]r ̅(x,y)=(r_0 ) ̅(x,y)+u ̅(x,y)[/math]. Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos [math]u ̅(ρ,ϑ)=〖(ρcosϑ-1)〗^2 ρ〖sin〗^2 ϑg ̅_ρ-〖(ρcosϑ-1)〗^2 sinϑcosϑg ̅_ϑ[/math]

2 Mallado

Comenzaremos planteando en Matlab el mallado de la placa con el fin de poder visualizarla en el espacio con mayor facilidad.

% Paso de muestreo 
h=0.1;
%Valores de la parametrización
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Conjunto de puntos que forman la malla
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujo del mallado
mesh(X,Y,0*X) 
title('Mallado')
%ajustar los ejes
axis([-3,3,-1,3])
view(2)

Figura 1

3 Temperatura

A continuación, escribiremos en matlab el codigo necesario para obtener la funcion temperatura de la placa T(x,y). Para ayudarnos a la visualización, colorearemos la placa con las distintas temperaturas en cada punto de la misma junto con las curvas de nivel, que unen los putnos con misma temperatura en la placa. Ademas, obtendremos la temperatura maxima en la placa. La funcion temperatura que definira el campo será: [math]T(x,y)=log⁡((y+x)^2+2)[/math] Nota: esta parte del codigo no es un programa en si. Para que funcione debidamente, debemos agregar la parte escrita en el apartado anterior.

%Campo de temperatura
 T=log(((Y+X).^2)+2);
 
 %Representación en 2D
subplot(1,2,1)
contour(X,Y,T)
view(2) 
axis([-3,3,-1,3]);
title('Curvas de nivel de la Temperatura')
%Representacion Temperatura 
subplot(1,2,1);
mesh(X,Y,T);
title('Temperatura')

subplot(1,2,2);
contour(X,Y,T,50);
colorbar
axis([-3,3,-1,3])
view(2)
title('Curvas de nivel')
%Temperatura máxima
Tmax=max(max(T));

4 Gradiente de temperatura

El gradiente de la temperatura nos definirá en que direccion se produce la maxima variacion de temperatura en la placa. Para poder mantener el orden dentro de la programacion, abriremos otro programa que repetirá exactamente los mismos comandos que en el apartado 1 de manera que se nos genere un mallado. En este caso le añadiremos el gradiente de la temperatura a toda la placa en forma de campo vectorial con flechas que nos indicará la direccion que sigue el gradiente. Estas flechas serán perpendiculares a las lineas de nivel que habiamos dibujado previamente en el apartado anterior.

%inicializamos las variables
 h=0.1;
 
 rho=[1:h:2];
 H=round(pi/h);         
 theta=linspace(0,pi,H);
 
 %Creamos el mallado y la parametrización
 [RHO,THETA]=meshgrid(rho,theta); 
 
 X=RHO.*cos(THETA);        
 Y=RHO.*sin(THETA);
 
 %Dibujamos el mallado
 mesh(X,Y,0*X) ; 
 axis([-3,3,-1,3]); 
 
 view(2);
 
 %GRADIENTE
 dx=((2*(Y+X))./(((Y+X).^2)+2));
 dy=dx;
 %CAMPO ESCALAR DE TEMPERATURAS
 T=log(((Y+X).^2)+2);
 %Creación de la subventana primera: gradiente en el campo de temperaturas
 
 subplot(1,2,1);
 hold on;
 surf(X,Y,T);
 
 %Dibujo del gradiente
 quiver3(X,Y,T,dx,dy,0*T);
 
 hold off;
 
 %Creación de la segunda subventana: campo vectorial del gradiente
 
 subplot(1,2,2);
 quiver(X,Y,dx,dy); 
 colorbar;

4.1 Variacion de temperatura

Calcularemos ahora la tasa de variacion de la temperatura que se produce al desplazarnos sobre la placa a una velocidad de 2m/s en la direccion del eje [math]i [/math]. Para esto, utilizaremos el gradiente calculado anteriormente y lo aplicaremos a la direccion antes mencionada, multiplicandola por la velocidad de desplazamiento dada. De esta forma, la tasa de variacion para este caso específico será calculada del siguiente modo:

1. Calculamos el gradiente de la temperatura(obtenido anteriormente): [math]∇T=(2*(Y+X)/((Y+X)^2)+2)[/math]
2. El vector velocidad sera el siguiente: [math]v=2i [/math]
3. Por ultimo aplicamos el gradiente a la direccion establecida para obtener la tasa: [math]TasaVariacion=∇T v =4/3[/math]

5 Campo de desplazamientos

Tomando el campo de desplazamientos [math]u ̅(ρ,ϑ)=〖(ρcosϑ-1)〗^2 ρ〖sin〗^2 ϑg ̅_ρ-〖(ρcosϑ-1)〗^2 sinϑcosϑg ̅_ϑ[/math], lo representaremos igual que el gradiente (flechas) con maor tamaño cuanta mayor es la intensidad del campo en una región.

%Paso de muestreo
h=0.1;
%Valores de parametrización 
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Conjunto de puntos que forman la malla
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 
view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
 %Componente x(X,Y) del campo vectorial
FX=zeros(31,16); 

 %Componente y(X,Y) del campo vectorial
FY=RHO.*(((cos(THETA)-1)).^2).*(RHO.*((sin(THETA)).^3))+(RHO.*sin(THETA).*cos(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY) 
axis([-3,3,-1,3])
view(2)
xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')

5.1 Campo de desplazamientos sobre la placa

Aplicaremos el campo de desplazamientos sobre la placa para poder apreciar la deformacion que sufre esta por una fuerza externa (generadora del campo de desplazamientos).

%Paso de muestro
h=0.1;
% Valores de parametrización
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Creación del mallado
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
%Componente x(X,Y) del campo vectorial
FX=zeros(31,16); 

 %Componente y(X,Y) del campo vectorial
FY=RHO.*(((cos(THETA)-1)).^2).*(RHO.*((sin(THETA)).^3))+(RHO.*sin(THETA).*cos(THETA));

%Matriz con las coordenadas en x, después del desplazamiento
Rx=FX+X;
%Matriz con las coordenadas en y, después del desplazamiento
Ry=FY+Y;

%Representación del sólido antes del desplazamiento 
subplot(1,2,1);
mesh(X,Y,0*X);
axis([-3,3,-1,8]);
title('Anillo sin deformación');
xlabel('Eje x');
ylabel('Eje y');

%Representación del sólido tras del desplazamiento
subplot(1,2,2);
mesh(Rx,Ry,0*Rx);
axis([-3,3,-1,8]);
title('Anillo con deformación');
xlabel('Eje x');
ylabel('Eje y');

6 Divergencia del campo de desplazamientos

Para obtener una medida de la variacion producida en la placa por el campo de desplazamientos, utilizaremos la divergencia de este. Aparecera en el grafico un mapa de colores con la variacion de volumen producida durante la deformación. Además, el programa nos devolverá el maximo valor que alcanza dicha divergencia.

%Paso de muestreo
h=0.1;
%Valores de parametrización 
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Conjunto de puntos que forman la malla
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 
view(2);
 
%Divergencia
DIVERGENCIA= [2.*((sin(THETA)).^2).*(RHO.*(cos(THETA)-1)).*(1+(RHO.*cos(THETA))+0.5.*(RHO.*cos(THETA)-1)+cos(THETA))]-[((cos(THETA)).^2).*((RHO.*cos(THETA))-1).^2];

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA) 
view(2)
colorbar
title('Divergencia')

7 Rotacional del campo de desplazamientos

Igual que con la divergencia, visualizaremos el rotacional del campo de desplazamientos con la siguiente programacion.

%Paso de muestreo
h=0.1;
%Valores de parametrización 
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Conjunto de puntos que forman la malla
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 
view(2);
 
%Rotacional
ROT= (-2.*(cos(THETA)-1./RHO).*sin(THETA).*cos(THETA).^2) - (2.*sin(THETA).*cos(THETA).*(RHO.*cos(THETA-1)).^2) + (2.*sin(THETA).^3) - (RHO.*(RHO.*cos(THETA)-1));

%Comprobar que está bien sacado el rotacional
%Representación del rotacional
surf(X,Y,ROT) 
view(2)
colorbar
title('Rotacional')

8 Tensiones

Definamos [math]e(u ̅ )=(∇u ⃗+〖∇u ⃗〗^t)/2[/math] , la parte simétrica del tensor gradiente de [math] u ̅[/math] conocido como tensor de deformaciones. En un medio elástico ideal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones a través de la fórmula [math] σ=λ∇∙u ̅1+2µe[/math], donde λ y µ son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. A pesar de que los desplazamientos son planos, (es decir, [math]u ̅ [/math]no tiene componente en la dirección [math]k ̅[/math]) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa. Tomando λ = µ = 1, dibujar las tensiones normales en la dirección que marca el eje [math]( g_ρ ) ⃗ [/math], es decir [math] (g_ρ ) ⃗∙σ∙(g_ρ ) ⃗[/math], las tensiones normales en la dirección que marca el eje [math] (g_θ ) ⃗/ρ[/math], es decir,[math] (g_θ ) ⃗/ρ∙σ∙(g_θ ) ⃗/ρ[/math], y las correspondientes al eje[math] k ⃗[/math], es decir,[math] k ⃗∙σ∙k ⃗[/math].

%Inicializamos las variables   
%Paso de muestreo
h=0.1;
%Valores de parametrización 
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Conjunto de puntos que forman la malla
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 
view(2);

%Gradiente de un campo vectorial
U=[(RHO.*cos(THETA)-1).^2.*RHO.*sin(THETA).^2 , -(RHO.*cos(THETA)-1).^2.*sin(THETA).*cos(THETA)];
[dx,dy]=gradient(U);

 %Dibujo de las tensiones respecto a g?
subplot(2,3,1);

 %Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
 r=dx(1,1);
 
 surf(X,Y,r)
 
 axis equal
 view(2)
 title('Tensiones respecto a g?')
 colorbar
 
 
 %Dibujo de las tensiones respecto a g?
 subplot(2,3,2);   
     
 %Tomamos el valor (2,2) de la matriz sigma y lo dibujamos   
 t=dy(2,2);
 
 surf(X,Y,t)
 
 axis equal
 title('Tensiones respecto a g?')
 view(2)
 
 colorbar
 
 %Dibujo de las tensiones respecto a gz
% subplot(2,3,3)
 
 %Tomamos el valor (3,3) de la matriz sigma y lo dibujamos 
 %z=           
     
 %surf(X,Y,z)
 
 %axis equal
% view(2)
% title('Tensiones respecto a gz')
 %colorbar

9 Tensión de Von Mises

La tensión de Von Mises se define con la formula [math]σ_VM=√((〖(σ_1-σ_2)〗^2+〖(σ_2-σ_3)〗^2+〖(σ_3-σ_1)〗^2)/2)[/math] donde [math]σ_1[/math],[math]σ_2[/math] y[math]σ_3[/math] son los autovalores de [math] σ [/math](también conocidos como tensiones superficiales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro). Pintar la tensión de Von Mises y señalar en que punto alcanza mayor valor. (Usar comando eig.m para calcular autovalores)

%Paso de muestreo
h=0.1;
%Valores de parametrización 
u=0.5:h:2; H=round(pi/h); v=linspace(0,pi,H);
%Conjunto de puntos que forman la malla
[RHO,THETA]=meshgrid(u,v); 
%Parametrización 
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);
 
 %Creamos la matriz sigma
 Sig=zeros(3,3);
 
 %Creamos la matriz de Von Mises 
 MatrizVonMises=zeros(16,21); 
 

 %Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
         for i=1:16
             for j=1:21
                r=RHO(i,j)';
                t=THETA(i,j)';
 
                Sig(1,1)=2*t/5+t*(2-1/r)/5;
 
                Sig(1,2)=(r-1)/(5*r^2);
 
                Sig(2,1)= (r-1)/(5*r^2);
 
                Sig(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
                Sig(3,3)=t*(2-1/r)/5;
 
              %Obtención de los autovalores de la matriz sigma
                [v,d]=eig(Sig); 
         
              %Fórmula de Von Mises
                MatrizVonMises11(i,j)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);  
 
            end
         end
 
     
 
 %Cambiamos a coordenadas cartesianas
 X=RHO.*cos(THETA);
 Y=RHO.*sin(THETA);
 
 %Dibujamos la tensión de Von Mises
 surf(X,Y,MatrizVonMises)
 
 colorbar
 view(2)
 title('Tensión de Von Mises')
 
 %Valor Máximo
 ValorMaximo=max(max(MatrizVonMises))
 
 %Coordenadas del valor Máximo
 [i,j]=find(MatrizVonMises==ValorMaximo))))
 
 RHO(i,j)
 THETA(i,j)

10 Masa de la placa

Finalmente, calcularemos la masa de la placa objeto de estudio. Esto se realizará mediante la integracion de la densida de la placa [math]d(x,y)=3+log⁡(1+|x|/(y+1)^4 )[/math]

%Inicializamos las variables    
h=0.1;
H=round(pi/h);

rho=[1:h:3];          
theta=linspace(0,pi,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);

d=3+log(1+abs(X)./(Y+1).^4);

% Representamos graficamente la densidad en cada punto de la malla
surf(X,Y,d);
axis([-3,3,-1,3]); 
axis square
%Cálculo de la masa total
MasaTotal=sum(sum(Densidad))