Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

Trabajo realizado por estudiantes
Título	Visualizaciones de campos escalares y vectoriales. Grupo 10-A
Asignatura	Teoría de Campos
Curso	2016-17
Autores	Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa, a la que se aplicará una temperatura T(x,y)=(y+2)·x², en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas. Por comodidad han sido elegidas las coordenadas cilíndricas. Los cálculo dispuestos estan formulados con matlab.

1 Mallado

Representaremos una malla de los puntos que componen la placa sólida con forma de corona de un cuarto de círculo. El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.

derecha

%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];          
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla 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([-1,4,-1,4]); 

view(2);
title('Mallado')

2 Temperatura

En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x²". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como curvas de nivel. El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

derecha

%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];          
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla 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 equal;

view(2);

%Campo escalar de temepraturas
Temperatura=(Y+2).*X.^2
TempMaxima=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==TempMaxima
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de las lineas de nivel en 2D 
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D 
subplot(1,2,2);

surf(X,Y,Temperatura);


view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')

3 Gradiente de la temperatura

Para conocer la variación de temperatura por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima. Se aprecia en la figura que el vector gradiente en perpendicular a la superficie de sólido.

derecha

%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

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

%Creamos el cojunto de puntos que forman la malla 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([-1,4,-1,4]); 
axis square

view(2);

%Gradiente
GradX=2*X.*Y+4*X;
GradY=X.^2;

%Campo escalar de temperaturas
Temperatura=(Y+2).*X.^2;

%Gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]); 
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]); 
axis square

hold off

%Campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY); 
axis([-1,4,-1,4]); 
axis equal
colorbar

4 Campo de desplazamientos

En este apartado procederemos a dibujar el campo de desplazamientos "u(ρ,θ)=θf(ρ)gρ". Las condiciones que conocemos son: los puntos de ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5, lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5. Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

derecha

%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

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

%Creamos el cojunto de puntos que forman la malla 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([-1,4,-1,4]); 

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA)); 

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')

5 Movimiento del sólido

En este caso se muestran dos gráficas una con la placa en reposo y otra después de aplicarle el campo de desplazamientos a la superficie, de esta forma es perfectamente apreciable la variación que sufre. El desplazamiento viene dado por dos componentes, el vector desplazamiento y los puntos de la placa respectivos a cada coordenada. Cuando ambas se suman, dan lugar a la figura final.

derecha

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

rho=[1:h:3];          
theta=linspace(0,pi/2,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([-1,4,-1,4]); 

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA)); 

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);  

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
% Movimiento de la placa en el eje X
MovimientoX =X+FX;

% Movimiento de la placa en el eje Y
MovimientoY =Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,X.*0)

axis([-1,4,-1,4]); 
view(2);
title('Aplicado el campo de movimientos')

6 Divergencia

Dada la divergencia usada para definir el campo "∇·(u)=θ(2-1/ρ)/5", procederemos a evaluar sus máximos, mínimos y nulos con el gradiente∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz. Al tener un dominio de ρ contenido en (1,3) y de θ en (0,pi/2), la derivada direccional gρ y gθ son positivas, también podemos sacar en claro que cuanto mayor sean ρ y θ mayor será su valor, llegando a su punto máximo en (3,pi/2) con: ∇·u(ρ,θ)=0.5236. A su vez el mínimo se encontrará donde los θ sean menores en nuestro dominio θ=0 dando lugar a: ∇·(u)=0 por estar ∇·(u) en función de θ.

derecha

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

rho=[1:h:3];          
theta=linspace(0,pi/2,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([-1,4,-1,4]); 

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')

7 Rotacional

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,theta)] = (rho-1)/(5*rho)gz. Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con theta=PI/4 y su valor es 0.1333.

El rotacional del campo de desplazamientos se ha calculado y aplicado a todos los puntos de la placa y han sido representados con colores.

derecha

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

rho=[1:h:3];          
theta=linspace(0,pi/2,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([-1,4,-1,4]); 

view(2);

%ROTACIONAL
ROT=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROT)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROT))

title('Rotacional')

8 Tensor de tensiones

Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

centro

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

centro

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

centro

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi.

centro

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

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

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

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

%En 2D

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

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5; 

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=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^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=THETA.*(2-1./RHO)./5;           
    
surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar



%Cálculo en 3D:
 

%Dibujo de las tensiones respecto a gρ 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gρ 3D')


%Dibujo de las tensiones respecto a gσ 3D
subplot(2,3,5);   
    
%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos   

surf(X,Y,t)

axis equal
title('Tensiones respecto a gσ 3D')


%Dibujo de las tensiones respecto a gz 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos            
    
surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')

centro

9 Tensiones tangenciales

Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica. Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)

10 Tensiones tangenciales respecto al plano ortogonal a ḡρ

derecha

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

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

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

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


%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gρ en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gρ en 2D')

11 Tensión de Von Mises

La tensión de Von Mises se define mediante la expresión:

centro

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 0.6283 alcanzándose e (1,pi/2).

derecha

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

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

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

%Creamos la matriz sigma
Sigma=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)';

               Sigma(1,1)=2*t/5+t*(2-1/r)/5;

               Sigma(1,2)=(r-1)/(5*r^2);

               Sigma(2,1)= (r-1)/(5*r^2);

               Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
               Sigma(3,3)=t*(2-1/r)/5;

             %Obtención de los autovalores de la matriz sigma
               [v,d]=eig(Sigma); 
        
             %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)

12 Masa total

La función de la densidad es: d(x,y)= 1 + e^-|x|/(y+1)2.

Una vez calculada, la masa total es m = 574,10 unidades

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

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

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

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

Densidad=1+exp(-abs(X)./(Y+1).^2);

% Representamos graficamente la densidad en cada punto de la malla
surf(X,Y,Densidad);
axis([-1,4,-1,4]); 
axis square

%Cálculo de la masa total
MasaTotal=sum(sum(Densidad))