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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualizaciones de campos escalares y vectoriales. Grupo 22-C
Asignatura Teoría de Campos
Curso 2015-16
Autores Vicente Elipe, Enrique Pírez, Angel Robisco, Juan Alfonso Oliveros, Santiago Valdés, Javier Olalde.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


Consideramos una placa plana que ocupa el anillo circular centrado en el origen y comprendido entre los radios 1 y 2 con el semiplano y>=0, en la cual definimos la temperatura \(T(x,y)\), dependiente de las variables \((x,y)\), y el campo de desplazamientos \(\vec u\) producido por la acción de una fuerza. De ellos, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) y el estudio del movimiento de la placa con el campo \(\vec u(\rho,\theta)\), con las consiguientes tensiones que producirá.

1 Mallado

Representaremos el mallado de los puntos que componen la placa sólida cuya forma es la de una corona semicircular. El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.

derecha
%Inicializamos las variables
h=0.1

rho=[1:h:2];          
theta=[0:h:pi];

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

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

%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-4,4,-4,4]); 

view(2);
title('Mallado')


2 Temperatura

La temperatura de la corona semicircular cambia según la función T(x,y)=log(y+2): La función de únicamente depende de la variable Y y se mantendrá constante para todo valor de X. Se ha dibujado el campo de temperaturas tanto en 2D como en 3D. El punto que ha resultado tener mayor temperatura T(x,y)=1,3863 es el punto x=0,y=2.

derecha
%Inicializamos las variables
h=0.1

rho=[1:h:2];          
theta=[0:h:pi];

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

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

%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-4,4,-4,4]); 

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=log(Y+2)
Maximo=max(max(Temperatura))

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

surf(X,Y,Temperatura);

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

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

surf(X,Y,Temperatura);
colorbar;
title('Temperatura en 3D')


3 Gradiente de la temperatura

El gradiente del campo de temperaturas (∇T) muestra la variación de la temperatura al movernos por el plano, siendo el módulo del gradiente lo que 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: donde la derivada direccional será máxima. El vector gradiente es perpendicular a la superficie del sólido. Cuanto más aumenta la Y, mayor es la temperatura: En el dibujo podemos observar el aumento de la temperatura a medida que los puntos del sólido avanzan en el eje Y.

derecha
%Inicializamos las variables
h=0.1

rho=[1:h:2];          
theta=[0:h:pi];

%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([-4,4,-4,4]); 

view(2);

%GRADIENTE
Grad=1./(Y+2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=log(Y+2);

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);

%Dibujo del gradiente
quiver3(X,Y,Temperatura,0*X,Grad,0*Temperatura);

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,0*X,Grad); 
colorbar


4 Campo de desplazamientos

El campo de desplazamientos dado es: u(rho,thetha)=1/5·Sin(2(thetha))*g(thetha). Vamos a dibujar el campo de vectorial, sabiendo que: g(thetha) = -rho·sen(thetha)i + rho·cos(thetha)j, cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY.

derecha
%Inicializamos las variables
h=0.1

rho=[1:h:2];          
theta=[0:h:pi];

%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([-4,4,-4,4]); 

view(2);

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

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

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

axis([-3,3,-3,3])
view(2)

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


5 Movimiento del sólido

Se ha representado el sólido antes y después del desplazamiento dado por el campo vectorial para así poder compararlos. Definiremos las funciones de desplazamiento en función de ambas componentes, sumando las componentes del vector de desplazamientos a las coordenadas respectivas de posición (Las de la posición incial). Se puede ver que el desplazamiento es mínimo.

derecha
%Inicializamos las variables
h=0.1

rho=[1:h:2];          
theta=[0:h:pi];

%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([-4,4,-4,4]); 

view(2);

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

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

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

surf(X,Y,0*X);  

axis([-4,4,-4,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,Y.*0)

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


6 Divergencia

Se ha calculado la divergencia del campo vectorial u en todos los puntos del sólido y la hemos representamos.

derecha
%Inicializamos las variables
rho=[1:0.1:2];          
theta=[0:0.1*pi:pi];

%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([-4,4,-4,4]); 

view(2);

%DIVERGENCIA
DIVERGENCIA=2/5*cos(2*THETA);

%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,thetha)] = (2/5)*Rho*Cos(2(Thetha)). Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos cuyo ángulo es pi/4.

derecha
%Inicializamos las variables
rho=[1:0.1:2];          
tetha=[0:0.1:pi];

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

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

%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-4,4,-4,4]); 

view(2);

%ROTACIONAL
ROTACIONAL=2/5*sin(2*THETA)

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

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

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

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)


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: 

                                           \(σ(ρ,θ)=λ∇.u.1 + 2με\

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

- Para hallar la dirección de ḡρ: (1,0,0)*σ*(1,0,0) = 2/5*cos(2θ)

- Para hallar la dirección de ḡθ: (0,1,0)*σ*(0,1,0) = 2/5*cos(2θ)*(1+(2/ρ^2))

- Para hallar la direcciñon de ḡz: (0,0,1)*σ*(0,0,1) = 2/5*cos(2θ)


Cálculo en 2D: 

%Inicializamos las variables   
h=0.1;

rho=1:h:2;
tetha=0:h:pi;

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

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

%Dibujo de las tensiones respecto a g sub rho
subplot(1,3,1);

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

surf(X,Y,r)

axis equal
view(2)
title('Dibujo de las tensiones respecto a g sub rho')
colorbar


%Dibujo de las tensiones respecto a g sub tetha
subplot(1,3,2);   
    
%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos   
t=(2/5.*cos(2*THETA)).*((1+2./RHO.^2));

surf(X,Y,t)

axis equal
title('Dibujo de las tensiones respecto a g sub tetha')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(1,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos 
z=2/5.*cos(2*THETA);           
    
surf(X,Y,z)

axis equal
view(2)
title('Dibujo de las tensiones respecto a g sub z')
colorbar


centro


Cálculo en 3D: 

%Inicializamos las variables   
h=0.1;

rho=1:h:2;
tetha=0:h:pi;

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

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

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(1,3,1);

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

surf(X,Y,r)

axis equal
title('Dibujo de las tensiones respecto a g sub rho 3D')
colorbar


%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(1,3,2);   
    
%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos   
t=(2/5.*cos(2*THETA)).*((1+2./RHO.^2));

surf(X,Y,t)

axis equal
title('Dibujo de las tensiones respecto a g sub tetha 3D')

colorbar

%Dibujo de las tensiones respecto a g sub z 3D
subplot(1,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos 
z=2/5.*cos(2*THETA);           
    
surf(X,Y,z)

axis equal
title('Dibujo de las tensiones respecto a g sub z 3D')
colorbar


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)

9.1 Tensiones tangenciales respecto al plano ortogonal a ḡρ

derecha
%Inicializamos las variables   
h=0.1;

rho=1:h:2;
tetha=0:h:pi;

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

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


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

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

surf(X,Y,Tang)

colorbar
title('Dibujo de las tensiones tangenciales respecto a g sub rho en 3D')

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

surf(X,Y,Tang)

colorbar
view(2)
title('Dibujo de las tensiones tangenciales respecto a g sub rho en 2D')


9.2 Tensiones tangenciales respecto al plano ortogonal a ḡθ

derecha
%Inicializamos las variables      
h=0.1;

rho=1:h:2;
theta=0:h:pi;

%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 theta
Tang=1./RHO.*(1/5.*sin(2*THETA)-1);

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

surf(X,Y,Tang)

colorbar
title('Dibujo de las tensiones tangenciales respecto a g sub theta en 3D')

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

surf(X,Y,Tang)

colorbar
view(2)
title('Dibujo de las tensiones tangenciales respecto a g sub theta en 2D')


10 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 2,0783.

derecha
%Inicializamos las variables    
h=0.1;

rho=[1:h:2];
theta=[0:h:pi];

%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(32,11); 


%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
    for i=1:32*11

        r=RHO(i)';

        t=THETA(i)';

        Sigma(1,1)=(2/5).*cos(2.*t);

        Sigma(1,2)=(1./r).*(0.2.*sin(2*t)-1);

        Sigma(2,1)=Sigma(1,2);

        Sigma(2,2)=Sigma(1,1)+(4/5).*cos((2*t)./(r.^2));

        Sigma(3,3)=Sigma(1,1);

        %Obtención de los autovalores de la matriz sigma
        [v,d]=eig(Sigma); 
        
        %Fórmula de Von Mises
        MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);  

    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
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 2,0783.


11 Masa total

La función de densidad es: d(x,y,z)=log(x^2+y^2+2). Hemos calculamos la masa y la hemos representado para visualizarlo de una manera más clara. Para hallar el valor de la integral hemos realizado la transformación de la función densidad a cilíndricas y hecho una aproximación utilizando el método del trapecio: la proyección del sólido se ha hecho tomando lado h=Particiones=1/1000 esperando lograr la mayor precisión posible.

La integral quedaría así:

centro

Una vez calculada, la masa total es m = 6.9975 

%Definimos el parámetro Rho entre los valores que toma
rho1=1
rho2=2

%Definimos la función de densidad en cilíndricas (Hay que acordarse del jacobiano de la transformación)
f=inline('rho.*log(rho.^2+2)')

%Definimos el número de particiones de la función
Particiones=1000

%Cálculo de la altura de los trapecios
Altura=(rho2-rho1)/Particiones; 

%Determinamos el número de valores que se van a evaluar
x=linspace(rho1,rho2,Particiones+1);

%Inicializamos el área
AREA=0;

    for i=1:Particiones

    %Determinamos el área de los trapecios diferenciales evaluando la función.
    A=(feval(f,x(i))+feval(f,x(i+1)))*Altura/2; 
    
    AREA=AREA+A;

    end
    
%El cálculo es una integral doble, pero como la función sólo depende de rho, se puede sacar la theta fuera.
Masa=pi.*AREA
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad._Grupo_22-C.&oldid=34043»