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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Deformaciones de una placa plana rectangular. Grupo C-1
Asignatura Teoría de Campos
Curso 2017-18
Autores Luis Pablo Rincón Bagüés, Adriana Mata Calvo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura



Consideramos una placa rectangular plana (en dimensión 2) que ocupa la region \((x,y) ∈ [-1,1]×[0,2]\). En ella vamos a suponer que tenemos definidas dos cantidades físicas:

  • La temperatura \(T(x,y)\), que depende de las variables espaciales \((x,y)\)
  • Los desplazamientos [math]\vec u(x,y) [/math] producidos por la acción de una fuerza determinada.

De esta forma, si definimos [math]\vec r_{0}(x,y) [/math] 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:

[math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y) [/math]

Vamos a suponer que la fuerza aplicada sonre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:

[math]\vec u(x,y) = x( -y + \frac{y^2}{2}) \vec i [/math].

1 Mallado

Dibujar un mallado que represente los puntos interiores del sólido. Tomar los ejes en el rectángulo \((x,y) ∈ [-2,2]×[-1,3] \) y como paso de muestreo [math] h=\frac{1}{10} [/math] para las variables \(x\) e \(y\).

Mallado del sólido
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(xx,yy);
%Dibujo 
mesh(xx,yy,0*xx)
%Ajuste de la gráfica
axis([-2,2,-1,3])
grid off
view(2)


2 Visualización campo de la temperatura

La temperatura del sólido viene dada por la función [math]T(x,y)=(x+2)(y-1)^2 +2[/math] en grados centígrados. Dibujar las curvas de nivel y decidir en qué punto la temperatura es máxima.
La función anterior es un campo escalar, es decir, es el campo de la temperatura del sólido.

Campo de la Temperatura
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
% TEMPERATURA (campo escalar)
T= inline('(x+2).*(y-1).^2+2','x','y');
zz=T(xx,yy);
% Dibujo
figure
 pcolor(xx,yy,zz);
 axis equal
 shading flat
 hold on
 contour(xx,yy,zz,'k')
 colorbar
 title('Temperatura máxima en (1,0)')
 hold off
% Ajuste de la gráfica
axis([-1,1,0,2])
view(2)


Se observa que en el punto \((1,0)\) la temperatura es máxima.

3 Cálculo y visualización del gradiente de Temperatura

El gradiente del campo de la temperatura (campo escalar) es el campo vectorial:

[math]\nabla T =\frac{\partial T}{\partial x_i} \vec e_i = (y-1)^2 \vec i + (x+2)(2y-2) \vec j [/math]

Como se puede observar en el dibujo el gradiente es ortogonal a las curvas de nivel del campo \(T\) y por tanto el [math]\nabla T[/math] es ortogonal a la superficie del campo \(T\) en todos los puntos de éste.

Campo vectorial: Gradiente de Temperatura
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
%Campo de Temperatura
T= inline('(x+2).*(y-1).^2+2','x','y');
%Gradiente
Tx= inline('(y-1).^2','x','y');
Ty= inline('(y.*2-2).*(x+2)','x','y');
%
zz=T(xx,yy);
gradTx=Tx(xx,yy);
gradTy=Ty(xx,yy);
% Dibujo de las curvas de nivel del campo de temperatura
figure 
quiver(xx,yy,gradTx,gradTy) 
title('Gradiente')
axis tight 
hold on 
Mz=T(xx,yy);
contour(xx,yy,zz,'r')
title('Curvas de nivel')
hold off


4 Variación de la temperatura

Calcular la variacion de temperatura por segundo si, partiendo del punto de coordenadas \((x,y)=(0,1)\) me muevo en la dirección \((1,1)\) con una velocidad de 2 metros por segundo.

  • Planteamiento:
  1. Partiendo del origen \(O\) nos movemos hacia \(P=(0,1)\), es decir, en la dirección [math] \vec e = \frac{\vec {OP}}{|\vec {OP}|} = \frac{\vec i}{|\vec i|} =\vec i [/math] con velocidad \(v=2\)m/s.
  2. De las propiedades del gradiente de un campo escalar:
[math] \frac{dT}{dt} = \frac{T(\vec r (t))}{dt} = \nabla{T(\vec r (t))} · \vec r ' = \nabla{T} · \vec v = \nabla{T} · (v \vec e) = v(\nabla{T} · \vec e)[/math]
y puesto que ya conocemos el \(\nabla{T}\):
[math] \nabla{T}(x,y) = (y-1)^2 \vec i + (x+2)(2y-2) \vec j [/math]
  • Entonces operando se obtiene:
[math] \frac{dT}{dt} = v·(y-1)^2 = 2(y-1)^2 [/math] ºC/s que en el punto \((1,1)\) es nula.

5 Visualización del campo de desplazamientos

A partir del vector de desplazamientos dado en el enunciado obtenemos el dibujo del campo de desplazamientos sobre los puntos del mallado del sólido.

Campo vectorial: Campo de desplazamientos
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
%Componentes en la dirección de i y de j del campo de desplazamientos 
ux=xx.*(-yy+(yy.^2)/2);
uy=0.*xx;
figure
%dibujo del mallado 
mesh(xx,yy,0*xx)
hold on 
%campo de desplazamientos
quiver(xx,yy,ux,uy,'k')
axis([-2,2,-1,3])
view(0,90)
hold off


6 Visualización del sólido antes y después del desplazamiento

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
%posicion final 
rx=(xx.*(-yy+yy.^2/2))+xx;
ry=0*xx+yy; 
%representacion de la superficie antes del desplazamiento  
figure
subplot(1,3,1)
surf(xx,yy,0*xx)
title('antes del desplazamiento')
axis([-2,2,-1,3])
view(2) 
xlabel('x')
ylabel('y') 
zlabel('z')
%representacion de la superficie después del desplazamiento
subplot(1,3,2) 
surf(rx,ry,0*rx)
axis([-2,2,-1,3])
title('después del desplazamiento')
axis([-2,2,-1,3])
xlabel('x')
ylabel('y') 
zlabel('z')
 
%comparacion
subplot(1,3,3)
surf(xx,yy,0*xx)
hold on 
surf(rx,ry,0*rx)
axis([-2,2,-1,3])
xlabel('x')
ylabel('y') 
zlabel('z') 
title('comparacion')
view(0,90)
hold off


7 Visualización divergencia del campo de desplazamientos

Representar el campo: [math]\nabla\cdot\vec{u}=-y+\frac{y^2}{2}[/math]

Visualización Campo Divergencia
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
%valor de la divergencia
divu=-yy+(yy.^2)/2;
%representacion 
surf(xx,yy,divu)
shading interp 
colorbar 
title('divergencia') 
xlabel('x')
ylabel('y') 
zlabel('z')


%Dibujo de la superficie que representa la divergencia y de sus curvas de
%nivel visto cuando z tiende a infinito 
figure
pcolor(xx,yy,divu)
shading interp 
hold on 
contour(xx,yy,divu,'k')
xlabel('x')
ylabel('y') 
hold off 
%Dibujo de la proyección sobre el plano YOZ de la superficie que representa
%la divergencia 
figure
plot(yy,divu)
title('proyección sobre el eje YOZ')
xlabel('y')
ylabel('z') 
grid on
maxi=max(max(divu)); 
mini=min(min(divu));


Visualización del campo escalar: Divergencia
Proyección de la Divergencia sobre el plano XOZ

No se produce aumento de volumen alguno, hay reducción de la parte central de la placa y los laterales no sufren alteraciones.
Como consecuencia de que la expresión de la [math]\nabla\cdot\vec{u}[/math] dependa únicamente de \(y\), ha sido posible determinar los puntos en los que la divergencia es máxima, mínima y nula gráficamente mediante la proyección de dicho campo sobre el plano \(YOZ\).

8 Rotacional

8.1 Cálculo del Rotacional

Siendo el campo de desplazamientos:

[math]\vec u[/math] como: [math]\vec u = u_x \vec i + u_y \vec j + u_z \vec k[/math]

ell rotacional de dicho campo se calcula como:

[math] \nabla\times \mathbf{u} =\left( \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z}\right)\vec i + \left(\frac{\partial u_x}{\partial z} - \frac{\partial u_z}{\partial x}\right)\vec j + \left(\frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y}\right)\vec k [/math]

Resultando finalmente:

[math]\nabla\times \mathbf{u}=x(1-y)\vec k [/math]

Siendo su módulo:

[math] |\nabla\times \mathbf{u}|= x(1-y)[/math]

8.2 Visualización

Visualización del rotacional
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
rot=xx.*(1-yy); 
%representacion del rotacional a partir del mallado 
surf(xx,yy,rot)
colorbar
axis([-2,2,-1,3])
view(2) 
title('rotacional') 
xlabel('x')
ylabel('y') 
zlabel('z')
%valores maximos y minimos 
rotmax=max(max(rot));
rotmin=min(min(rot));

Los puntos que sufren un mayor rotacional son \((x,y)=(-1,2)\) y \((x,y)=(1,0)\) y el menor son \((x,y)=(-1,0)\) y \((x,y)=(-2,1)\).

9 Tensiones Normales

9.1 Cálculo del tensor de deformaciones

Se llama tensor de deformaciones a la parte simétrica del tensor gradiente de [math]\vec u [/math]

[math]e(\vec u) =\frac{(\nabla{\vec u}+\nabla{\vec u}^t)}{2}[/math]

Para ello debemos calcular [math]\nabla{\vec u}[/math] y [math]\nabla{\vec u}^t[/math]

  • [math]\nabla{\vec u}=\frac{\partial{\vec u}}{\partial{x_m}}\otimes\vec e_m =(-y+\frac{y^2}{2})(\vec i\otimes\vec i)+(x(-1+y))(\vec i\otimes\vec j) [/math]
  • [math]\nabla{\vec u}^t =(-y+\frac{y^2}{2})(\vec i\otimes\vec i)+(x(-1+y))(\vec j\otimes\vec i)[/math]

Finalmente obtenemos el tensor de deformaciones:

  • [math]e(\vec u)=(-y+\frac{y^2}{2})(\vec i\otimes\vec i)+\frac{x}{2}(-1+y)[(\vec i\otimes\vec j)+(\vec j\otimes\vec i)][/math]

9.2 Cálculo del tensor de tensiones [math]\sigma_{ij}[/math]

En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{ij}[/math] a través de la fórmula

[math]\sigma = λ\nabla · \vec u 1 +2μe [/math]

donde λ y μ son 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]\vec u[/math] no tiene componente en la dirección de [math]\vec k[/math]) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa.

Tras haber calculado el tensor de deformaciones y la divergencia del campo [math]\vec u[/math] en los apartados anteriores, y tomando λ=μ=1 calculamos el tensor de tensiones

[math]\sigma = (-y+\frac{y^2}{2}) 1 + 2[(-y+\frac{y^2}{2})(\vec i\otimes\vec i)+\frac{x}{2}(-1+y)[(\vec i\otimes\vec j)+(\vec j\otimes\vec i)][/math]

9.3 Cálculo de las tensiones normales

9.3.1 Tensiones normales en la dirección que marca el eje [math]\vec i [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec{i}[/math] se obtienen a partir de la expresión:

[math]\vec{i}\cdot\sigma_{ij}\cdot\vec{i}=3(-y+\frac{y^2}{2})[/math]

9.3.2 Tensiones normales en la dirección que marca el eje [math]\vec j [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec{j}[/math] se obtienen a partir de la expresión:

[math]\vec{j}\cdot\sigma_{ij}\cdot\vec{j}=(-y+\frac{y^2}{2})[/math]

9.3.3 Tensiones normales en la dirección que marca el eje [math]\vec k [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec{k}[/math] se obtienen a partir de la expresión:

[math]\vec{k}\cdot\sigma_{ij}\cdot\vec{k}=(-y+\frac{y^2}{2})[/math]

9.4 Visualización de las tensiones normales

Tensiones normales 

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
 
%tensiones normales en las direcciones de los ejes coordenados sobre el
%mallado 
TNi=(-yy+(yy.^2)/2).*3;
TNj=-yy+(yy.^2)/2;
TNk=-yy+(yy.^2)/2; 
%representacion de las tensiones normales 
figure
subplot(1,3,1) 
surf(xx,yy,TNi)
axis image
view(2) 
colorbar
xlabel('x')
ylabel('y') 
zlabel('z')
title('Tensiones normales en la dirección del eje i')
 
subplot(1,3,2) 
surf(xx,yy,TNj)
axis image
view(2) 
colorbar
xlabel('x')
ylabel('y') 
zlabel('z')
title('Tensiones normales en la dirección del eje j')
 
subplot(1,3,3)
surf(xx,yy,TNk)
axis image
view(2) 
colorbar
xlabel('x')
ylabel('y') 
zlabel('z')
title('Tensiones normales en la dirección del eje k')

10 La Tensión de Von Mises

La tensión de Von Mises se define por la fórmula:

[math]\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}[/math]

Donde [math]\sigma_{1}, \sigma_{2}[/math] y [math]\sigma_{3}[/math] son los autovalores de [math]\sigma[/math] (también conocidos como tensiones principales) cuya expresión es:

[math](\sigma_{ij})=\begin{pmatrix} -3y+\frac{3y^2}{2} & -x+xy & 0 \\ -x+xy & -y+\frac{y^2}{2} & 0 \\ 0 & 0 & -y+\frac{y^2}{2} \end{pmatrix}[/math]

La tensión de Von Mises es una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro). 

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
%calculo de los autovalores 
A= inline('3.*(-y+(y.^2)./2)','x','y');
B= inline('x.*(-1+y)','x','y');
C= inline('-y+(y.^2)./2','x','y');
sigma=[];
 for i=1:length(x)
     for j=1:length(y)
         J=A(xx(i,j),yy(i,j));
         K=B(xx(i,j),yy(i,j));
         I=C(xx(i,j),yy(i,j));
         sigma=[J K 0; K I 0 ; 0 0 I ];
         AU=eig(sigma);  %autovalores de sigma 
         VM=sqrt(((AU(1)-AU(2))^2+(AU(2)-AU(3))^2+(AU(3)-AU(1))^2)/2);
         zz(i,j)=VM;
     end 
 end 
%dibujo 
figure
subplot(1,2,1)
surf(xx,yy,zz)
colorbar
xlabel('x')
ylabel('y') 
zlabel('z')
title('Tensión de Von Mises')
%Estudio del máximo valor de la tensión de Von Mises
subplot(1,2,2)
hold on
m=zz(1,1:length(zz));
%Dibujo de la proyección sobre el plano XOZ de la superficie que genera la tensión de Von Mises
plot(x,m,'b');
maxzz=max(m);
%Bucle que pinta los puntos donde la tensión de Von Mises es máxima sobre la proyección
for k=1:length(m);
   if m(k)==maxzz
         plot(x(k),maxzz,'xr','markersize',10)
     end
end
axis([-3 3 -0 2]);
xlabel('Eje x')
ylabel('Eje z')
legend('Proyección sobre el plano XOZ','Punto máximo')
hold off
fprintf('El máximo valor de la tensión de Von Mises es %f\n',maxzz)

Tensión de Von Mises y proyección sobre el plano XOZ El máximo valor de la tensión de Von Mises se alcanza para [math]x=±1[/math] y el valor de la tensión de Von Mises para dicho punto es de \(1.732051\).



11 Campo de Fuerzas

El campo de fuerzas [math]\vec F[/math] que actúa sobre la placa (y que son causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal

[math]\vec{F}=-\nabla\cdot\sigma=-\frac{\partial \sigma_{ji}}{\partial x_{j}}\vec{e}_{i}[/math]

Lo calculamos

  • [math]\vec F_1=-\frac{\partial \sigma_{j1}}{\partial x_{j}} \vec e_{1} = -[\frac{\partial (-y+\frac{3y^2}{2})}{\partial x} + \frac{\partial(-x+xy)}{\partial y} + \frac{\partial 0}{\partial z}]\vec i = -x \vec i[/math]
  • [math]\vec F_2=-\frac{\partial \sigma_{j2}}{\partial x_{j}} \vec e_{2} = -[\frac{\partial (-x+xy)}{\partial x} + \frac{\partial (-y+\frac{y^2}{2})}{\partial y} + \frac{\partial 0}{\partial z}]\vec j = 2(1-y) \vec j [/math]
  • [math]\vec F_3=-\frac{\partial \sigma_{j3}}{\partial x_{j}} \vec e_{3} = -[\frac{\partial 0}{\partial x} + \frac{\partial 0}{\partial y} + \frac{\partial (-y+\frac{y^2}{2})}{\partial z}]\vec k = 0 \vec k[/math]

Entonces, la expresión del campo de fuerzas [math]\vec{F}[/math] es:

[math]\vec{F}= -x \vec i + 2(1-y) \vec j [/math]
Campo de fuerzas
%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
%componentes del campo de fuerzas
Fx=-xx; 
Fy=2.*(1-yy);
%Dibujo del campo de fuerzas 
figure
mesh(xx,yy,0*xx) 
hold on 
quiver(xx,yy,Fx,Fy,'k')
axis image
view(2) 
xlabel('x')
ylabel('y')
title('campo de fuerzas')
hold off


  • Interpretación de la gráfica

En el dibujo podemos ver un campo de fuerzas cuyo modulo decrece progresivamente a medida que nos acercamos al centro de esta simetría radial en el punto \((0,1)\) en donde se anula.

12 Cálculo de la masa total de la placa

Suponiendo que la placa tiene una densidad que responde a la ecuación:

[math]d(x,y)=3+e^{\frac{-|x|}{(y+1)^{4}}}[/math]

Calculamos la masa total de la placa:

[math]\int_{A} d(x,y) dA=\int_{-1}^{1}\int_{0}^{2} d(x,y) dydx=\int_{-1}^{1}\int_{0}^{2} 3+e^{\frac{-|x|}{(y+1)^{4}}} dydx=15.71 [/math]

Sin embargo, puesto que, en esta ocasión, la primitiva es difícil de calcular, se recurre a la utilización de métodos numéricos. Por esta razón, empleamos la fórmula del trapecio para la aproximación de integrales: 

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-1:0.1:1;
y=0:0.1:2;
%Creación del mallado
[xx,yy]=meshgrid(x,y);  
%extremos integracion: 
a=-1; b=1; c=0; d=2; 
N1=length(x)-1; N2=N1; %subintervalos 
h1=(b-a)/N1; h2=(d-c)/N2; 
u=a:h1:b; 
v=c:h2:d; 
[U,V]=meshgrid(u,v);
%densidad 
densidad=3+exp(-abs(U)./((V+1).^4));
%formula del trapecio 
w1=ones(N1+1,1);
w1(1)=1/2; 
w1(N1+1)=1/2;
w2=ones(N2+1,1);
w2(1)=1/2;
w2(N2+1)=1/2; 
masa=h1*h2*(w2')*densidad*w1;
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad._Grupo_C-1&oldid=41468»