Visualización de campos escalares y vectoriales en elasticidad en placa plana

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en elasticidad en placa plana. Grupo 2-C
Asignatura Teoría de Campos
Curso 2019-20
Autores Joaquín del Campo Lobo, Sergio Martínez Alcántara, Javier Montero Barrionuevo, Daniel Cavallé Pulla
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 región (x,y) ∈ [−2,2]×[0,4]. En ella supondremos que tenemos definidas dos cantidades físicas: la Temperatura T(x,y) como campo escalar, y los desplazamientos [math]\vec u(x,y)[/math] como campo vectorial, 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 vendrá dada por [math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)[/math] Supondremos que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos [math]\vec u(x,y) = xf(y) \vec j [/math] donde f(y) es una cierta función.

1 Representación del mallado del sólido

Definimos la placa mediante dos vectores x e y que representan el intervalo en el cual está definida nuestra placa, tomando un muestreo h=1/10. Al crear la malla donde la placa está definida, al representarla con el comando surf, tenemos que tener en cuenta que es plana dando valor 0 a la altura de los puntos de la misma. 

h=0.1;                       %muestreo
x=-2:h:2;                    %eje x de la placa
y=0:h:4;                     %eje y de la placa
[xx,yy]=meshgrid(x,y);       %mallado de la placa
figure(1)
surf(xx,yy,0*xx);            %representación de la placa
axis([-3,3,0,5])             %ejes del rectángulo
view(2)
Mallado de la placa

2 Representación del campo escalar. Temperatura

El campo escalar viene dado por la temperatura [math]T(\rho,\theta)=\rho^2sin(\theta)[/math] que viene dada en coordenadas polares. Realizando el cambio a cartesianas para su mejor visualización y operabilidad dada la geometría de la placa, [math] \rho =\sqrt{x^2 + y^2} [/math] y [math] \rho sin(\theta) =y [/math] , nos queda la temperatura como [math]T(x,y)=y\sqrt{x^2 + y^2} [/math]. Aplicando el campo a nuestra placa nos queda la siguiente visualización:

Aplicación del campo T en la placa

obtenida con el código: 

h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
f=sqrt(xx.^2+yy.^2).*yy
surf(xx,yy,f)
axis([-3,3,0,5])
view(2)
colorbar


Para la mejor visualización del campo de Temperatura en la placa, dibujamos las líneas de nivel del mismo, obteniendo el máximo valor en el punto (0,4) de la placa con un valor de 17.889.

Curvas de nivel del campo T en la placa

obtenida con el código: 

h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
f=inline('yy.*sqrt(xx.^2+yy.^2)','xx','yy')
zz=f(xx,yy)
[C,h]=contour(xx,yy,zz,14)
axis([-3,3,0,5])
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2);
 title('FUNCIÓN DE TEMPERATURA');                       
xlabel('EJE X');                                        
ylabel('EJE Y');                                        
zlabel('EJE Z');                                        
M=max(max(T))                                             
colormap hot;
colorbar ;


Ahora pasamos a calcular el gradiente del campo de Temperaturas ([math]\nabla T[/math])

Para el cálculo de [math]\nabla T[/math] en coordenadas cartesianas: [math]\nabla T=\frac { \partial T(x,y) }{ \partial x } \overrightarrow { i } \quad +\frac { \partial T(x,y) }{ \partial y } \overrightarrow { j } \quad [/math].

Para nuestro campo de temperaturas [math]T(x,y)=y\sqrt{x^2 + y^2} [/math], obtenemos:
[math]\nabla T=\frac {xy}{\sqrt{x^2 + y^2}}\overrightarrow { i } \quad + \frac {y^2}{\sqrt{x^2 + y^2}}\overrightarrow { j } \quad+\sqrt{x^2 + y^2} \overrightarrow { j } \quad \longrightarrow [/math] [math] \nabla T=\frac {xy}{\sqrt{x^2 + y^2}} \overrightarrow { i } \quad + \frac {2y^2+x^2}{\sqrt{x^2 + y^2}}\overrightarrow { i } \quad [/math]

Como podemos observar en el gráfico, las curvas de nivel del campo son ortogonales a la dirección del gradiente.

gradiente
h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
fx=xx.*yy./sqrt(xx.^2+yy.^2);
fy=2*yy.^2+xx.^2/sqrt(xx.^2+yy.^2);
f=inline('yy.*sqrt(xx.^2+yy.^2)','xx','yy');
T=f(xx,yy);
quiver(xx,yy,fx,fy);
hold on
contour(xx,yy,T)
axis([-3,3,0,5])
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


3 Representación del campo vectorial. Desplazamientos

Consideramos ahora el campo de desplazamientos [math]\vec u [/math] del cual sabemos que:

  • Los puntos situados en el eje y=0 no sufren desplazamiento en la dirección [math]\vec j [/math].
  • El [math]\nabla·\vec u =\frac{x}{10}[/math].

Por consiguiente el campo es del tipo [math]\vec u =xf(y)\overrightarrow { j } [/math]

Aplicando la definición de la divergencia en coordenadas cartesianas:

[math]\nabla·\vec u(x,y) =\frac{\partial}{\partial x^i}({u_i})=\frac{\partial}{\partial y}({xf(y)})=\frac{x}{10} \longrightarrow [/math] [math] xf'(y) = \frac{x}{10} \longrightarrow [/math] [math] f'(y)= \frac{1}{10}\longrightarrow [/math] [math] f(y)=\frac{y}{10}+C[/math]

Con la primera de las condiciones obtenemos el valor de la constante, C=0, ya que [math] f(0)=0[/math]

Por lo que el campo de desplazamientos es: [math]\vec u = \frac{xy}{10}\overrightarrow { j }[/math]

Dibujamos el el campo de vectores en los puntos del mallado del sólido con el código: 

h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
mesh(xx,yy,0*xx)
axis([-3,3,0,5])
figure(1)
fx=yy.*0./10;
fy=xx.*yy./10;
quiver(xx,yy,fx,fy);
view(2)


campo de vectores

4 Deformación de la placa

Como podemos observar en las gráficas, la placa sufre los mayores desplazamientos en los laterales de la misma, siendo el borde izquierdo (x=-2) el que se comprime, mientras que el borde derecho (x=2) se alarga. Estas deformaciones coinciden con el campo de deformaciones [math] \vec u [/math] que hemos dibujado en el apartado anterior.

sólido antes y después del desplazamiento


sólido antes y después del desplazamiento


h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
axis([-3,3,0,5]);
subplot(1,2,1)                        
mesh(xx,yy,0*xx)
axis([-3,3,0,5]);
title('ANTES DEL DESPLAZAMIENTO');    
xlabel('EJE X');                      
ylabel('EJE Y');                      
subplot(1,2,2) 
fx=yy.*0./10;
fy=xx.*yy./10;           
mesh(xx+fx,yy+fy,0*xx)
axis([-3,3,0,5]);
title('DESPUÉS DEL DESPLAZAMIENTO'); 
xlabel('EJE X');                     
ylabel('EJE Y');


5 Interpretación de la divergencia y del rotacional

Como podemos observar en la gráfica de la divergencia, la placa sufrirá, como nos dice la definición de divergencia, una variación de volumen local debido al campo de desplazamientos. Dicha variación, como bien se puede apreciar en la gráfica, es máxima en los extremos de la placa, de sentido opuesto. Puesto que la divergencia es una función lineal [math]\nabla·\vec u =\frac{x}{10}[/math], tanto el máximo y el mínimo de la función se encuentran en los extremos del intervalo en el cual está definida la función ([-2,2] en nuestro caso). La divergencia es nula donde se anula la divergencia, en x=0.

gráfica de la divergencia


Pasamos ahora a interpretar el rotacional del campo [math]\vec u [/math]:

[math]\nabla \times \vec u=\left|\begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ & & \\ u_x & u_y & u_z \end{matrix}\right| =\left| \begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ & & \\ 0 & \frac{xy}{10} & 0 \end{matrix}\right| =\left( - \frac{ \partial \frac{xy}{10}}{\partial z}\right)\vec i + \left(\frac{ \partial \frac{xy}{10}}{\partial x} \right)\vec k = \frac{y}{10}\vec k [/math]

Por lo tanto su módulo será igual a: [math]|\nabla×\vec u| =\frac{y}{10}[/math]. Como apreciamos en la gráfica, los puntos que sufren un mayor rotacional son aquellos que se encuentran en el borde de la placa (y=4) ya que el rotacional al depender sólo de y, es máximo en el extremo y=4.

grafica del rotacional
h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
mesh(xx,yy,0*xx)
axis([-3,3,0,5]);
view(2)
figure(1)
f=(yy./10)
surf(xx,yy,f)
axis([-3,3,0,5]);
view(2)
colorbar


6 Tensiones normales

Definimos [math]\epsilon (\overrightarrow { u } )=\frac { \nabla \overrightarrow { u } +{ \nabla }\overrightarrow { u } ^{ t } }{ 2 } [/math], conocido como tensor de deformaciones. Para calcular este tensor, que se define como la parte simétrica del tensor gradiente de [math]\ \vec u, \nabla \vec u[/math], debemos sacar el gradiente de [math]\ \vec u[/math] y su traspuesto [math]\ \vec u^t[/math]. Puesto que nuestro campo de desplazamientos es [math]\vec u = \frac{xy}{10}\overrightarrow { j }[/math]:

[math]\nabla \overrightarrow { u } =\left( \begin{matrix} \frac { \partial { \overrightarrow { u } }_{ 1 }}{\partial x} & \frac { \partial { \overrightarrow {u} }_{ 1 } }{ \partial y } & \frac { \partial { \overrightarrow { u } }_{ 1 } }{ \partial z } \\ \frac { \partial { \overrightarrow { u } }_{ 2 } }{ \partial x } & \frac { { \partial \overrightarrow { u } }_{ 2 } }{ \partial y } & \frac { { \partial }\overrightarrow { u } _{ 2 } }{ \partial z } \\ \frac { \partial { \overrightarrow { u } }_{ 3 } }{ \partial x } & \frac { { \partial \overrightarrow { u } }_{ 3 } }{ \partial y } & \frac { { \partial \overrightarrow { u } }_{ 3 } }{ \partial z } \end{matrix} \right) [/math].

obtenemos [math]\nabla \vec u = \left( \begin{matrix} 0 & 0 & 0 \\ \frac { y }{ 10 } & \frac {x}{10} & 0 \\ 0 & 0 & 0 \end{matrix} \right)[/math] y [math]\nabla \vec u^t = \left( \begin{matrix} 0 & \frac { y }{ 10 } & 0 \\ 0 & \frac {x}{10} & 0 \\ 0 & 0 & 0 \end{matrix} \right)[/math]

Por lo que el tensor de deformaciones [math]\epsilon ({ \overrightarrow { u } }) =\left( \begin{matrix} 0 & \frac { y }{ 20 } & 0 \\ \frac { y}{ 20 } & \frac {x}{10} & 0 \\ 0 & 0 & 0 \end{matrix} \right) [/math].


Una vez hallado el tensor de deformaciones, se puede hallar el tensor de tensiones mediante la expresión [math]\sigma_{ij} = \lambda \nabla \cdot \vec u \boldsymbol {1} + 2 \mu \epsilon [/math]. La divergencia [math]\ \nabla · \vec u[/math] ya la tenemos calculada: [math] \nabla·\vec u =\frac {x}{10}[/math]. Sean [math] \lambda[/math] y [math] \mu [/math] los Coeficientes de Lamé, ambos los consideramos de valor 1, utilizando el cálculo tensorial obtenemos [math]\sigma_{ij}[/math]:


[math] \sigma_{ij} = \lambda \nabla · \vec u \boldsymbol {1} + 2 \mu \epsilon = \frac {x}{10} ( \vec i \otimes \vec i + \vec j \otimes \vec j + \vec k \otimes \vec k ) + 2 ( \frac {y}{20} \vec i \otimes \vec j + \frac {y}{20}\vec j \otimes \vec i + \frac {x}{10}\vec j \otimes \vec j) = \frac {x}{10} \vec i \otimes \vec i + \frac {3x}{10} \vec j \otimes \vec j + \frac {x}{10} \vec k \otimes \vec k +\frac {y}{20} \vec i \otimes \vec j + \frac {y}{20}\vec j \otimes \vec i[/math]


[math] \sigma_{ij} = \begin{pmatrix} \frac {x}{10} & \frac {y}{10} & 0\\ \frac {y}{10} & \frac {3x}{10} & 0\\ 0 & 0 & \frac {x}{10} \end{pmatrix} [/math]

Para calcular las tensiones normales en las direcciones de los ejes, basta con tomar los elementos de la diagonal:
[math] \vec i · \sigma_{ij} · \vec i = \frac {x}{10} [/math] en la dirección que marca [math] \vec i[/math].
[math] \vec j · \sigma_{ij} · \vec j = \frac {3x}{10} [/math] en la dirección que marca [math] \vec j[/math].
[math] \vec k · \sigma_{ij} · \vec k = \frac {x}{10} [/math] en la dirección que marca [math] \vec k[/math].

Dibujamos las tensiones normales:

tensiones normales en la dirección i
h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
mesh(xx,yy,0*xx)
sigmai=((xx).*(1/10))
mesh(xx,yy,sigmai)
axis([-3,3,0,5]);
title('TENSIONES NORMALES EJE i (2D)')   
xlabel('EJE X')                          
ylabel('EJE Y')          
colorbar;


tensiones normales en la dirección j
h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
mesh(xx,yy,0*xx)
sigmai=((3.*xx).*(1/10))
mesh(xx,yy,sigmai)
axis([-3,3,0,5]);
title('TENSIONES NORMALES EJE j (2D)')   
xlabel('EJE X')                          
ylabel('EJE Y')          
colorbar;


tensiones normales en la dirección k
h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
mesh(xx,yy,0*xx)
sigmai=((xx).*(1/10))
mesh(xx,yy,sigmai)
axis([-3,3,0,5]);
title('TENSIONES NORMALES EJE k (2D)')   
xlabel('EJE X')                          
ylabel('EJE Y')          
colorbar;


Para calcular las tensiones tangenciales al plano ortogonal a [math]\ \vec i [/math], que viene dada por [math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math].
Operando, las tensiones tangenciales [math] \sigma_{ti} = \frac {y}{10}[/math].

7 Tensión de Von Misess

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 [/math], [math] \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores de [math] \sigma [/math] (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.

miniaturadeimagen


h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
mesh(xx,yy,0*xx);
fx=inline('xx./10','xx','yy');
fxy=inline('yy./10','xx','yy');
fxz=inline('0*xx','xx','yy');
fyx=inline('yy./10','xx','yy');
fy=inline('3.*xx./10','xx','yy');
fyz=inline('0.*xx','xx','yy');
fzx=inline('0.*xx','xx','yy');
fzy=inline('0.*xx','xx','yy');
fz=inline('xx./10','xx','yy');
for j= 1:length(x)                              % Inicio del bucle
for i= 1:length(y)                                                                              % Asignación de valores para cada componente del mallado                    
Fx=fx(xx(i,j),yy(i,j));                           % Respecto a fx
Fxy=fxy(xx(i,j),yy(i,j));
Fxz=fxz(xx(i,j),yy(i,j));
Fyx=fyx(xx(i,j),yy(i,j));
Fy=fy(xx(i,j),yy(i,j));                       % Respecto a fy
Fyz=fyz(xx(i,j),yy(i,j));
Fzx=fzx(xx(i,j),yy(i,j));
Fzy=fzy(xx(i,j),yy(i,j));
Fz=fz(xx(i,j),yy(i,j));               
sigma=[Fx,Fxy,Fxz;Fyx,Fy,Fyz;Fzx,Fzy,Fz]; 
autov=eig(sigma)
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(xx,yy,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 valor máximo es igual a 0.81888, que se alcanza en el punto [-2,4], esquina superior izquierda de la plaza.

8 Campo de fuerzas [math]\ \vec F [/math] sobre la placa

El desplazamiento de la placa provocado por el vector [math]\vec u(x, y) = \frac{xy}{10} \vec j [/math] está producido por un vector de fuerza [math]\ \vec F [/math] que actúa en el interior de la placa.

El cálculo de dicho vector de fuerzas se realiza haciendo una aproximación de la ecuación de la elasticidad lineal, que es : [math]\ \vec F = \nabla ·\sigma_{ij}= -\frac {\partial \sigma_{ij}}{\partial x_j} · \vec e_i[/math]
Obteniendo las componentes de dicho vector aplicando la divergencia a la matriz: [math]\ \sigma_{ij}[/math].

[math]\ \vec F_i = -\frac {\partial \sigma_{j1}}{\partial x_j} · \vec e_i = -\left( \frac {\partial \frac {x}{10}}{\partial x} + \frac {\partial \frac {y}{10}}{\partial y} + \frac {\partial 0}{\partial z} \right) · \vec i = \frac {-1}{5} \vec i[/math]

[math]\ \vec F_j = -\frac {\partial \sigma_{j2}}{\partial x_j} · \vec e_j = -\left( \frac {\partial \frac {y}{10}}{\partial x} + \frac {\partial \frac {3x}{10}}{\partial y} + \frac {\partial 0}{\partial z} \right) · \vec j = 0 \vec j[/math]

[math]\ \vec F_k = -\frac {\partial \sigma_{j3}}{\partial x_j} · \vec e_k = -\left( \frac {\partial 0}{\partial x} + \frac {\partial 0}{\partial y} + \frac {\partial \frac {x}{10}}{\partial z} \right) · \vec k = 0 \vec k[/math]

Dibujando el campo de fuerzas [math]\ \vec F [/math] podemos observar que solamente actúa en la dirección [math]\ \vec i [/math] (ya que es la única componente del campo). Esto provoca el campo de desplazamientos [math]\ \vec u [/math] en la dirección normal al campo de Fuerzas, que como hemos calculado antes, dicho campo de desplazamientos actúa en la dirección [math]\ \vec j [/math]

Campo de fuerzas F
h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
mesh(xx,yy,0*xx)
axis([-3,3,0,5])
figure(1)
fx=-xx./5;
fy=0;
quiver(xx,yy,fx,fy);
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
view(2)


9 Cálculo de la masa de la placa

Para calcular la masa de la placa nos viene indicada la densidad de la misma, que es [math] d(x,y) = \frac {1}{x^2 + y^2 +2} [/math] Utilizando la integración numérica para su cálculo, el valor de la masa de la placa es: [math] M = 30.855 [/math] (unidades másicas). 

h=0.1;
x=-2:h:2;
y=0:h:4;
[xx,yy]=meshgrid(x,y);
figure(1)
mesh(xx,yy,0*xx)
d=1+exp(-abs(xx)./(yy+1).^2);
mesh(xx,yy,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
dm=abs(d*h^2);                        % Masa de cada cuadrado que conforma el mallado
masa=sum(sum(dm))


densidad de la masa
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad_en_placa_plana&oldid=48876»