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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo 1-A
Asignatura	Teoría de Campos
Curso	2016-17
Autores	Samantha Valencia Zapata Belén Cisneros Araujo Inas Rahmouni Sefiani Alicia Jiménez Barral
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Consideramos una placa rectangular plana (dimensión 2) que ocupa la región . En ella vamos a suponer que tenemos definidas dos cantidades fisicas: la temperatura T(x,y), que depende de las dos variables espaciales (x,y), y los desplazamientos del vector producidos por la acción de una fuerza determinada. De esta forma, si definimos el vector de 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:

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 donde f(y) es una cierta función.

2 Mallado de la placa

Representamos la placa rectangular en un mallado con la función "mesh". El mallado está definido en la región . Al ser el sólido plano damos valor nulo a las cotas. La gráfica tiene representados unos intervalos en los ejes de coordenadas. Definimos estos intervalos con la función "axis". En el eje de abscisas el intervalo es [-3,3] y en el de ordenadas es [0,5] con un muestreo de altura 0,1.

x=-2:0.1:2;                                      % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                                       % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                             % Obtención del mallado 
mesh(X,Y,0*X)                                    % Representación superficial con líneas entrecruzadas
axis([-3,3,0,5])                                 % Establecimiento de valores mínimos y máximos en los ejes representados
title('MALLADO DE LA PLACA RECTANGULAR PLANA')   % Título del gráfico
xlabel('EJE X')                                  % Título del eje x
ylabel('EJE Y')                                  % Título del eje y

3 Temperatura de la placa

Definimos la temperatura del sólido con la función:

La temperatura es un campo escalar, es decir, que a cada punto (x,y) se le asocia un escalar T(x,y) de una región determinada. Este campo escalar se representa mediante las superficies de nivel, que son el lugar geométrico de los puntos del espacio para los cuales la función escalar toma un valor constante. Éstas al ser cortadas por un plano se convierten en curvas de nivel, denominadas isotermas.

Definimos la función temperatura con la función "inline" y las curvas de nivel con la función "contour". En la gráfica podemos observar que las isotermas aumentan de temperatura según se alejan de la parte inferior izquierda de la placa. Por consiguiente, el punto donde la temperatura es máxima es el (2,4), al que le corresponde la temperatura de valor 52.

x=-2:0.1:2;                                             % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y =0:0.1:4;                                             % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y] = meshgrid(x,y)                                   % Obtención del mallado
ft=inline('((x+2).^2+(y+2).^2)','x','y');               % Definimos la función temperatura.     
T = ft(X,Y) 
[C,h]=contour(X,Y,T,14);                                % Definimos las curvas de nivel.
axis([-3,3,0,5]);                                       % Establecimiento de valores mínimos y máximos en los ejes representados
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2);
title('FUNCIÓN DE TEMPERATURA');                        % Título del gráfico
xlabel('EJE X');                                        % Título del eje x
ylabel('EJE Y');                                        % Título del eje y
zlabel('EJE Z');                                        % Título del eje z
max(max(T))                                             % Obtención del valor máximo.
colormap hot;
colorbar ;

4 Gradiente de temperatura

Dada la función de temperatura T(x,y) anterior, se define su vector gradiente como las derivadas parciales respecto x e y de la función T(x,y). Por lo tanto, se calculan dichas derivadas parciales: y .

El campo vectorial del gradiente se representa con la función "quiver". A continuación, se superponen las curvas de nivel con "hold on" y "hold off".

El gradiente es perpendicular a las superficies de nivel. Esta propiedad queda reflejada en el gráfico, ya que los vectores inciden ortogonalmente en las curvas de nivel.

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)    
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
                                         % Definimos las componentes del campo del gradiente de la temperatura.
Tx=2*X+4;
Ty=2*Y+4;
                                         % Representamos el campo vectorial generado por el gradiente
quiver(X,Y,Tx,Ty)
hold on                                  % Superposición del gradiente con las lineas de nivel
contour(X,Y,T)                           % Representación de las lineas de nivel
axis([-3,3,0,5]);                        % Establecimiento de valores mínimos y máximos en los ejes representados
title('CAMPO GRADIENTE DE TEMPERATURA'); % Título del gráfico
xlabel('EJE X');                         % Título del eje x
ylabel('EJE Y');                         % Título del eje y 
hold off

5 Campo de desplazamiento

Consideramos el campo de desplazamientos u: Cuya divergencia es:

La divergencia de un campo vectorial es un campo escalar. Cuya fórmula en coordenadas cartesianas es:

Para definir la función f(y) se calcula la divergencia del vector de desplazamientos y se iguala a la divergencia que nos dan en el enunciado. Despejando f ’(y) obtenemos que es igual a 1/10. Después integramos f ’(y) respecto dy, obteniendo que f(y)=y/10 + cte. Como dato especifican que los puntos en el eje y=0 no sufren desplazamiento en la dirección j. Por lo tanto, sustituimos y=0 en la ecuación anterior y nos queda que la constante es nula.

x=-2:0.1:2;                                   % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                                    % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                          % Obtención del mallado
U=inline('(-X.^3)./20 + (Y./10)','x','y');    % Definimos el campo de vectores U y sus componentes
Ux=inline('(-x.^3)./20','x');
Uy=inline('y./10','y');
U=Ux(X);
V=Uy(Y);
quiver(X,Y,U,V)                               % Dibujamos el campo vectorial
axis([-3,3,0,5]);                             % Establecimiento de valores mínimos y máximos en los ejes representados
title('CAMPO DE DESPLAZAMIENTOS');            % Título del gráfico
xlabel('EJE X');                              % Título del eje x
ylabel('EJE Y');                              % Título del eje y

6 Campo de vectores de desplazamiento

El campo de desplazamientos finalmente es:

Representamos dicho campo sobre el mallado con la función "quiver", estableciendo anteriormente y .

Al observar el gráfico comprobamos que los vectores situados en el eje x no se mueven en la dirección j ya que y=0. Según se va incrementando el valor de y, los puntos empiezan a sufrir el desplazamiento en esa dirección.

x=-2:0.1:2;                                     % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                                      % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                            % Obtención del mallado
hold on                                         % Superposición del gradiente con las lineas de nivel
mesh(X,Y,0*X)                                   % Representación superficial con líneas entrecruzadas
axis([-3,3,0,5])                                % Establecimiento de valores mínimos y máximos en los ejes representados
                                                % Definición del punto de observación(azimut=0,elevación=90)
xlabel('EJE X')                                 % Título del eje x
ylabel('EJE Y')                                 % Título del eje y
title('CAMPO DE DESPLAZAMIENTOS EN EL MALLADO') % Título del gráfico 
Ux=(-X.^3)./20;
Uy=(Y./10);
quiver(X,Y,Ux,Uy)                               % Dibujamos el campo vectorial
hold off

7 Placa antes y después del desplazamiento

Los desplazamientos del vector están producidos por la acción de una fuerza determinada. En primer lugar representamos la placa antes de la deformación mediante el vector de posición donde podemos observar en la primera gráfica la forma que tiene inicialmente.

Dicha fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos: donde f(y) ha sido calculada anteriormente y toma el valor de y/10.

A continuación, sabiendo que cada punto (x,y) adquiere una nueva posición obtenemos la representación de una segunda placa. La posición de cada punto (x,y) de la placa después de la deformación viene dada por la suma del vector de posición inicial y el vector de desplazamientos:

Una vez representada la placa en los dos gráficos mediante el comando "subplot", se observa que ha sufrido deformación tanto en el sentido longitudinal como en el sentido transversal, produciéndose un alargamiento en el primer caso y un estrechamiento en el segundo.

x=-2:0.1:2;                           % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                            % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                  % Obtención del mallado
subplot(1,2,1)                        % Representación del sólido antes del desplazamiento
mesh(X,Y,0*X)
axis([-3,3,0,5])                      % Establecimiento de valores mínimos y máximos en los ejes representados
title('ANTES DEL DESPLAZAMIENTO');    % Título del gráfico
xlabel('EJE X');                      % Título del eje x
ylabel('EJE Y');                      % Título del eje y
subplot(1,2,2)                        % Representación del sólido después de desplazarse
Ux=-X.^3/20;
Uy=Y./10;
mesh(X+Ux,Y+Uy,0*X)

axis([-3,3,0,5]);                    % Establecimiento de valores mínimos y máximos en los ejes representados
title('DESPUÉS DEL DESPLAZAMIENTO'); % Título del gráfico
xlabel('EJE X');                     % Título del eje x
ylabel('EJE Y');                     % Título del eje y

8 Divergencia del desplazamiento

La divergencia es una medida del cambio de volumen local debido al desplazamiento que se produce en los puntos de la placa. Nos piden determinar analíticamente los puntos donde la divergencia de u es máxima, mínima y nula. Para ello vamos a obtener su representación en 2D y 3D mediante la función "mesh". La divergencia del campo de desplazamientos u viene dada por:

En las gráficas podemos observar la forma que adquiere la representación de la divergencia. Si la proyectamos en el plano XZ vemos que tiene la forma de una parábola donde su vértice es el punto donde la divergencia es máxima. Los puntos donde la divergencia es mínima se encuentran en la intersección de la parábola con el plano z=-0.5 y los puntos donde es nula se encuentran en la intersección de la parábola con el plano z=0.

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
div=(-3*X.^2)./20+1./10;                 % Definimos la divergencia
subplot(1,2,1)                           % Representación 2D    
mesh(x,y,div);                            
axis([-3,3,0,5])                         % Establecimiento de valores mínimos y máximos en los ejes representados 
title('DIVERGENCIA DESPLAZAMIENTOS')     % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y

colorbar;
subplot(1,2,2)                           % Representación 3D
mesh(x,y,div);
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y
zlabel('EJE Z')                          % Título del eje z
title('INTERPRETACIÓN ESPACIAL')         % Título del gráfico

colorbar;
axis([-3,3,0,5,-0.5,0.1])                % Establecimiento de valores mínimos y máximos en los ejes representados

9 Rotacional del desplazamiento

El rotacional sobre campos vectoriales muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto. El rotacional se calcula mediante la fórmula ,siendo F el vector desplazamiento, si lo calculamos, obtenemos un resultado nulo. Es coherente a la hora de ver la gráfica, ya que en ningún punto muestra la tendencia a inducir rotación.

10 Tensor de deformaciones

En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones σ_ij a través de la fórmula:

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 el vector u no tiene componente en la dirección de k) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa.

Definimos ɛ como la parte simétrica del tensor gradiente del vector u conocido como tensor de deformaciones:

Tomando λ=µ=1, se pide dibujar las tensiones normales en la dirección que marca el eje i, es decir i·σ·i, las tensiones normales en la dirección que marca el eje j, es decir j·σ·j y las correspondientes al eje k, es decir k·σ·k.

En primer lugar comenzamos definiendo ɛ para poder calcular el valor de σ. Para ello necesitamos calcular el gradiente de un campo vectorial mediante la siguiente expresión:

Una vez obtenido el valor del gradiente del vector u se suma a su traspuesta y se divide entre dos para obtener el valor de ɛ para sustituirlo en la fórmula del tensor de tensiones. Por lo que sustituyendo también la divergencia del vector u y tomando λ=µ=1, el valor de σ obtenido es el siguiente:

Ahora se procede a calcular las tensiones normales en la dirección que marcan los ejes i, j y k.

10.1 Tensiones normales en la dirección del eje i

La tensiones normales en la direccion del eje i vienen expresadas por:

Se representan en 2D y 3D:

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
sigmai=(-9.*(X.^2)./20)+ 1./10;          % Definición del campo de tensiones
subplot(1,2,1)                           % Representación 2D
mesh(X,Y,sigmai)
axis([-3,3,0,5])                         % Establecimiento de valores mínimos y máximos en los ejes representados
title('TENSIONES NORMALES EJE i (2D)')   % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y

subplot(1,2,2)                           % Representación 3D
mesh(X,Y,sigmai)
title('TENSIONES NORMALES EJE i (3D)')   % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y
zlabel('EJE Z')                          % Título del eje z

10.2 Tensiones normales en la dirección del eje j

La tensiones normales en la direccion del eje j vienen expresadas por:

Se representan en 2D y 3D:

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
sigmaj=(-3.*(X.^2)./20)+ 3./10;          % Definición del campo de tensiones
subplot(1,2,1)                           % Representación 2D
mesh(X,Y,sigmaj)
axis([-3,3,0,5])                         % Establecimiento de valores mínimos y máximos en los ejes representados
title('TENSIONES NORMALES EJE j (2D)')   % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y

subplot(1,2,2)                           % Representación 3D
mesh(X,Y,sigmaj)
title('TENSIONES NORMALES EJE j (3D)')   % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y
zlabel('EJE Z')                          % Título del eje z

10.3 Tensiones normales en la dirección del eje k

La tensiones normales en la dirección del eje k vienen expresadas por:

Se representan en 2D y 3D:

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
sigmak=(-3.*(X.^2)./20)+ 1./10;          % Definición del campo de tensiones
subplot(1,2,1)                           % Representación 2D
mesh(X,Y,sigmak)
axis([-3,3,0,5])                         % Establecimiento de valores mínimos y máximos en los ejes representados
title('TENSIONES NORMALES EJE k (2D)')   % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y

subplot(1,2,2)                           % Representación 3D
mesh(X,Y,sigmak)
title('TENSIONES NORMALES EJE k (3D)')   % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y
zlabel('EJE Z')                          % Título del eje z

11 Tensiones tangenciales

A partir del tensor de tensiones σ calculado en el apartado anterior, podemos hallar también las tensiones tangenciales respecto al plano ortogonal ῑ. Para ello utilizaremos la siguiente expresión:

Como podemos ver tras estos cálculos llegamos a la conclusión de que dichas tensiones son nulas.

12 Tensión de Von Mises

La tensión de Von Mises se define por la fórmula: donde σ1, σ2 y σ3 son los autovalores de σ (también conocidos como tensiones principales). 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). Nos pide pintar la tensión de Von Mises y señalar en qué punto se alcanza el mayor valor.

El código Matlab elaborado para la representación es el siguiente:

x=-2:0.1:2;                                     % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                                      % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                            % Obtención del mallado 
                                                % Definimos la función 
fx=inline('(-9*x.^2)./20+1/10','x','y');
fy=inline('(-3*x.^2)./20+3/10','x','y');
fz=inline('(-3*x.^2)./20+1/10','x','y');

for i= 1:length(x)                              % Inicio del bucle
  for j= 1:length(x)                             
                                                % Asignación de valores para cada componente del mallado                    
xx=fx(X(i,j),Y(i,j));                           % Respecto a fx
yy=fy(X(i,j),Y(i,j));                           % Respecto a fy
zz=fz(X(i,j),Y(i,j));                           % Respecto a fz
v=[xx yy zz];                                   % Creación de un vector con las componentes [xx, yy, zz]
N=diag(v);                                      % Diagonalización de ese vector
autov=eig(N);                                   % Obtención de los autovalores
                                                % Cálculo de la tensión de Von Mises para cada componente

valormaximo=sqrt(((autov(1)-autov(2))^2+((autov(2)-autov(3))^2+((autov(3)-autov(1))^2))*1/2));  
Z(i,j)=valormaximo;
     end
end                                             % Fin del bucle
axis([-3,3,0,5])                                % Establecimiento de valores mínimos y máximos en los ejes representados
mesh(X,Y,Z)                                     % Representación superficial con líneas entrecruzadas
title('TENSIÓN DE VON MISES')                   % Título del gráfico
xlabel('EJE X')                                 % Título del eje x
ylabel('EJE Y')                                 % Título del eje y
zlabel('EJE Z')                                 % Título del eje z
maximo=max(max(valormaximo))                    % Obtención del valor máximo

El máximo valor que alcanza es 1.5620

13 Campo de fuerzas

El campo de fuerzas F que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal: donde la matriz σ calculada anteriormente es:

Para representar el campo de fuerzas F aplicamos la ecuación para cada componente de la matriz σ. Como resultado la ecuación es nula en todas las componentes excepto en la primera cuyo valor es:

Esto significa que las fuerzas solo actúan en la dirección i. Si observamos la placa después de la deformación vemos que sufre un estrechamiento transversal que hace que adopte una forma mas alargada. Esto se debe a que las fuerzas que actúan en sentido contrario a los dos lados de la placa, la comprime, dando lugar a la deformación citada anteriormente.

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
Fx=inline('X.*(18./20)','X');            % Definimos el campo de fuerzas F y sus componentes
Fy=0.*Y;
U=Fx(X);
V=Fy;
quiver(X,Y,U,V)                          % Dibujamos el campo de fuerzas
title('CAMPO DE FUERZAS')                % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y
axis([-3,3,0,5])                         % Establecimiento de valores mínimos y máximos en los ejes representados

14 Masa de la placa

La función de densidad de la placa viene expresada por:

Primero vamos a representar dicha función para ver como se distribuye la densidad en la placa:

x=-2:0.1:2;                              % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:0.1:4;                               % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                     % Obtención del mallado
d=1+exp((-(abs(X)))./((Y+1).^2));        % Definimos la función densidad
mesh(X,Y,d);
title('FUNCIÓN DENSIDAD');               % Título del gráfico
xlabel('EJE X');                         % Título del eje x
ylabel('EJE Y');                         % Título del eje y
zlabel('EJE Z');                         % Título del eje z

Para calcular la masa de la placa a partir de la función de densidad hemos utilizado el método del Trapecio. Este consiste en calcular aproximadamente el valor de una integral definida dividiendo la proyección del sólido en cuadrados de lado h=1/1000 para conseguir una mayor precisión. La masa total de la placa es la suma del producto de la densidad de la placa y del área de cada diferencial de masa.

h=1/1000;                             % Lado de los cuadrados en los que dividimos el sólido
x=-2:h:2;                             % Vector x con valores entre -2 y 2 con saltos de h (1/10)
y=0:h:4;                              % Vector y con valores entre 0 y 4 con saltos de h (1/10)
[X,Y]=meshgrid(x,y);                  % Obtención del mallado
d=1+exp((-(abs(X)))./((Y+1).^2));     % Valor de la densidad para cada punto del mallado
dm=abs(d*h^2);                        % Masa de cada cuadrado que conforma el mallado
masa=sum(sum(dm))                     % Suma de las masas de todos los cuadrados

La masa total de la placa es: 29.4891