Visualización de campos escalares y vectoriales en elasticidad (Grupo A-4)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad (Grupo A-4)
Asignatura	Teoría de Campos
Curso	2016-17
Autores	María Guadalupe Aranda Sánchez, Ignacio Blázquez González, Javier Alejandro Chirivella Colosso, María Lucía Cueva Duarte, Javier Palomares Alonso
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Enunciado

En este trabajo vamos a representar y analizar campos escalares y vectoriales en el espacio en relación a su elasticidad.

Para ello consideramos una placa plana rectangular en dos dimensiones. La placa está acotada por x=-2 y x=2, y por y=0 e y=4.

Para representar la placa utilizaremos coordenadas cartesianas. Para ello hemos pasado los datos del enunciado que nos venían en coordenadas polares a cartesianas.

Tenemos

y

Entonces:

Por otro lado, consideramos el campo de desplazamiento ū, sabiendo:

    -ū=x·f(y)j
    -
    -Los puntos del eje y=0 no se desplazan en el sentido de j

Por otro lado nos dan la fórmula para hallar la tensión a partir del desplazamiento:

donde ε(ū)=(grad ū+grad ū^t)/2, y tomamos λ=μ=1 (coeficientes de Lamé)

También la tensión de Von Mises:

donde σ₁, σ₂ y σ₃ son los autovalores de σ.

y el campo de fuerzas F:

Finalmente,

sabemos que la densidad de la placa es d(x,y)= 1 + e^-|x|/(y+1)2.

2 Creación del mallado

En primer lugar dibujaremos un mallado que contenga todos los puntos del interior del sólido. Tomamos los ejes en el rectángulo definido por [math] -3≤x≤3;0≤y≤5 [/math], y utilizamos como paso de muestreo [math] h=1/10 [/math] para las variables x e y. Para ello utilizamos el siguiente código:

%1.Crear un mallado que represente los puntos del sólido

%Paso
h=0.1;
%Definimos los vectores posicion de la placa
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
Mz=0.*Mu;
%Creamos la placa
axis([-3,3,0,5])
figure(1)
surf(Mu,Mv,Mz)
view(2)

Obtenemos un mallado, que será la base sobre la que trabajaremos:

Mallado inicial

3 Representación de la temperatura

Seguidamente representamos, dentro de ese mallado, la función temperatura, que depende de las variables x e y. Como ya hemos dicho, la función temperatura es:

Para representarla, dibujamos las curvas de nivel de la función dentro del sólido. En nuestro caso, hemos decidido representar 100 curvas de nivel.

%Curvas de nivel Temperatura

%Paso
h=0.1;
%Definimos los vectores posicion de la placa
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
Mw=Mv.*(sqrt(Mu.^2+Mv.^2));
%Creamos la placa
axis([-3,3,0,5])
figure(1)
surf(Mu,Mv,Mw)
contour3(Mu,Mv,Mw,100)
view(2)

Curvas de nivel de la función temperatura

Para deducir dónde es máxima la temperatura, acudimos a la siguiente imagen en 3 dimensiones, generada por Matlab:

Representación 3D. Se ven con facilidad los máximos

Como podemos ver en la imagen, las temperaturas máximas se dan en los puntos (-2,4) y (2,4).

4 Cálculo del gradiente de temperaturas

A continuación, calculamos el gradiente de temperaturas y lo representamos como un campo vectorial.

%Paso
h=0.1;
%Definimos los vectores posicion de la placa
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
V2=((Mu.^2).*(Mv.^2).*Mv)./sqrt((Mu.^2)+(Mv.^2));
V1=Mu.*Mv./sqrt((Mu.^2)+(Mv.^2));
quiver(Mu,Mv,V1,V2)
view(2)

Obtenemos el campo vectorial del gradiente:

Gradiente de temperaturas

Podemos comprobar como el gradiente es ortogonal a las curvas de nivel superponiendo las dos representaciones:

Gradiente de temperaturas y Curvas de Nivel

5 Cálculo del campo de desplazamientos

Esta vez se nos pide obtener el campo de desplazamientos ū, con las condiciones citadas en el enunciado:

    -ū=x·f(y)j
    -
    -Los puntos del eje y=0 no se desplazan en el sentido de j

Con estas condiciones podemos calcular el campo ū analíticamente:

ū=x·y/10

6 Campo de vectores de los puntos del mallado

Ahora lo que haremos es representar el campo de vectores de los puntos del mallado calculado anteriormente. Para ello elaboramos el siguiente código:

%Paso
h=0.1;
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
%Definimos e campo de vetores u.
V1=0*Mu;
V2=Mu.*Mv./10;
quiver(Mu,Mv,V1,V2)
view(2)

Y el resultado que sale por pantalla es el siguiente:

Representación de los campos de los puntos del mallado

7 Sólido antes y después del desplazamiento

Seguidamente tenemos el sólido antes de la aplicación del desplazamiento, así como el resultado de aplicarlo después. Aplicamos el comando subplot, como se nos indica, para que las dos representaciones queden sobre la misma gráfica. El código es:

h=0.1;
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
%Definimos el campo de vectores u.
V1=0*Mu;
V2=Mv;
fx=0;
fy=(Mu.*Mv)./10;
Mz=0.*Mu;
hold on
axis([-3,3,0,5])
subplot(2,1,1);
surf(Mu,Mv,0*Mu) 
view(2)
subplot(2,1,2);
surf(Mu+fx,Mv+fy,0.*Mu)
view(2)
hold off

Y obtenemos el gráfico:

Representación del sólido antes y después del desplazamiento

8 Representación la Divergencia de ū

La divergencia nos sirve para analizar cambios locales de volumen producidos por el desplazamiento.

Para calcularlo utilizamos el siguiente código:

h=0.1;
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
Div=Mu/10;
surf(Mu,Mv,Div)

La imagen que obtenemos es:

Representación dela divergencia de ū

Podemos observar como la divergencia tiene relación con la dirección del campo de desplazamiento de cada uno de los puntos, de forma que es negativa hasta x=0,donde el volumen se reduce, y positiva a partir de x=0, ya que el volumen aumenta. En x=0, la divergencia también es 0 pues no hay desplazamiento.

9 Cálculo del módulo del rotacional de ū en los puntos del sólido

Calculamos el módulo del rotacional de ū utilizando el siguiente código:

%Paso
h=0.1;
u=[-2:h:2];
v=[0:h:4];
w=u.*0;
[Mu,Mv,Mw]=meshgrid(u,v,w);
%Definimos el campo de vectores u.
V1=0.*Mu;
V2=0.*Mu;
V3=Mv./10;
quiver3(Mu,Mv,Mw,V1,V2,V3)

Y obtenemos el siguiente gráfico:

Representación del módulo del rotacional en los puntos del sólido

Podemos observar que los valores máximos del rotacional se alcanzan cuando y=-4.

10 Representación las tensiones normales en la dirección del eje X, del eje Y y del eje Z

Se nos da la fórmula para hallar la tensión a partir del desplazamiento:

donde ε(ū)=(grad ū+grad ū^t)/2, y tomamos λ=μ=1 (coeficientes de Lamé).

Primero escribimos el código para calcular la tensión en el eje X:

h=0.1;                 
u=-2:h:2;               
v=0:h:4;           
[Mu,Mv]=meshgrid(u,v);
%TENSIONES NORMALES DEL EJE X
ty=Mv*0;
tx=Mu/10;
% dibujamos %  
quiver(Mu,Mv,tx,ty)
ylabel('tensiones normales eje x')

Y para el eje Y:

h=0.1;                 
u=-2:h:2;               
v=0:h:4;           
[Mu,Mv]=meshgrid(u,v);
%TENSIONES NORMALES DEL EJE Y
ty=3*Mu/10;
tx=0*Mv;
% dibujamos %  
quiver(Mu,Mv,tx,ty)
ylabel('tensiones normales eje y')

Para el eje Z no hay tensión, así que no lo consideramos.

Las representaciones resultantes son:

Representación de la Tensión en la dirección OX

Representación de la Tensión en la dirección OY

Si combinamos las dos representaciones en un solo gráfico obtenemos:

Representación de las Tensiones en las direcciones OX (azul) y OY (verde)

11 Cálculo de tensiones tangenciales

Calculamos las tensiones tangenciales respecto al plano ortogonal a i, es decir |σ·ī-(ī·σ·ī)ī|.

El código es:

h=0.1;                 
 u=-2:h:2;               
 v=0:h:4;           
 [Mu,Mv]=meshgrid(u,v);
%TENSIONES TANGENCIALES DEL EJE X
 ty=Mu*0;
 tx=Mv/10;
 % dibujamos %  
 quiver(Mu,Mv,tx,ty)
 ylabel('tensiones tangenciales eje x')

Y obtenemos:

Representación de las Tensiones tangentes al plano ortogonal a ī

12 Representación de la tensión de Von Mises

Utilizamos la fórmula:

donde σ₁, σ₂ y σ₃ son los autovalores de σ.

El código es el siguiente:

h=.1;
u=[-2:h:2];
v=[0:h:4];
[Mu,Mv]=meshgrid(u,v);
A=[];
for  i=1:length(u)
    for j=1:length(v)
    A=[Mu(i,j)./10 Mv(i,j)./10 0;Mv(i,j)./10 (3.*Mu(i,j))./10 0;0 0 Mu(i,j)./10];  
    L=eig(A);
    ecuacion= sqrt(((L(1)-L(2))^2+(L(2)-L(3))^2+(L(3)-L(1))^2)/2); 
    Mw(i,j)=ecuacion;
    end                                 
end
surf(Mu,Mv,Mw)

Obtenemos la representación:

Representación de las Tensiones de Von Mises

13 Dibujar el Campo de Fuerzas F

Se nos da la fórmula:

Se trata de una aproximación del campo de fuerzas que actúan sobre la placa, y que son las causantes del desplazamiento, mediante la ecuación de la elasticidad.

Utilizamos el código:

%F=[-1/5 0 0]
%Paso
h=0.1;
u=[-2:h:2];
v=[0:h:4];
w=u*0;
[Mu,Mv,Mw]=meshgrid(u,v,w);
%Definimos e campo de vetores u.
V1=-Mu/10-Mv/10;
V2=Mu*0;
V3=Mw*0;
quiver3(Mu,Mv,Mw,V1,V2,V3)

Y nos devuelve el siguiente gráfico:

Representación del campos de fuerzas F

14 Cálculo de la masa total de la placa

Para finalizar, calculamos la masa de la placa, sabiendo que la densidad de la placa es d(x,y)= 1 + e^-|x|/(y+1)2.

El código es el siguiente:

h=0.1;                 
u=-2:h:2;               
v=0:h:4;           
[Mu,Mv]=meshgrid(u,v);
Ds=(Mu.*Mv).*(1./(2+(Mu.^2)+(Mv.^2)));
ms=abs(Ds.*h^2);
Masa=sum(sum(ms))

Obtenemos