Visualización de campos escalares y vectoriales en elasticidad (Grupo B2)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad
Asignatura	Teoría de Campos
Curso	2022-23
Autores	David Astorga Trejos Abdelali Zariohi Boutaleb Sergio Míguez González Nacira Faraji Bahja Luis David Sánchez Pérez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa rectangular plana (en 2D) que ocupa la región [math](x, y) ∈ [0, 10] × [0, 2][/math]. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x, y)[/math], que viene dada por,

[math]T(x, y)=(x-6)^2+(10(y-\frac{3}{2})^2[/math]

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 [math](x, y)[/math] 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 sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos

[math]\vec{u}(x, y,t)=\vec{a} sin(\vec{d}*\vec{r}-vt))\vec{i},[/math]

donde [math]\vec{a}[/math] se conoce como amplitud, [math]k\gt0[/math],es el número de onda, [math]\vec{d}[/math] es un vector unitario que marca la dirección de propagación y v es la velocidad de propagación.

Para nuestro trabajo, tendremos los siguientes valores: [math] \vec{a}=\frac{2}{5}\vec{i}[/math], [math]\vec{d}=\vec{i}[/math] y [math]k=1 [/math]

Operando con estos valores obtenemos el siguiente vector [math]\vec{u}(x,y,t)[/math]

[math]\vec{u}(x,y,t)=\frac{2}{5} sin(x)\vec{i}[/math]

1 Mallado de la placa.

Primero debemos representar el sólido definido en el enunciado acorde a los parámetros especificados en el enunciado. Las directrices requeridas son las siguientes: "Tomar los ejes (comando axis) en el rectángulo [math](x, y) ∈ [−0.5; 10.5] × [−0.5; 2.5][/math] y como paso de muestreo [math]h = \frac{2}{10}[/math] para las variables x e y.

Representación del mallado.

%parametrización con muestreo h=2/10.
r=linspace(0,10,50);
t=linspace(0,2,10);

%mallado
[R,T]= meshgrid(r,t);      
          
Z=zeros(10,50);         
surf(R,T,Z)
%ejes
axis([-0.5,10.5,-0.5,2.5])
view(2)
%nombramos los ejes
xlabel('Eje X')
ylabel('Eje Y')

2 Temperatura de la placa.

En este apartado realizaremos un diagrama de temperaturas de la placa que represente por colores la variación de la temperatura y otro para dibujar las líneas de dicho campo escalar.

Representación de la temperatura y sus correspondientes curvas de nivel.

r=linspace(0,10,50);
t=linspace(0,2,10);
[x,y]= meshgrid(r,t);      
          
T=((x-6).^2)+(10.*(y-3/2)).^2;

hold on
%División de la ventana en dos
subplot(1,2,1)              
 %Representamos el campo escalar de temperaturas
surf(x,y,T)                  
view(2)
axis([-0.5,10.5,-0.5,2.5])
xlabel('Eje X')
ylabel('Eje Y')
colorbar                      %Mostramos la escala
subplot(1,2,2)                %Trabajamos en la segunda ventana
contour(x,y,T,30)
title('Curvas de Nivel','Fontsize',16)      %Líneas de nivel
colorbar                      %Mostramos las escala
axis([-0.5,10.5,-0.5,2.5])
xlabel('Eje X')
ylabel('Eje Y')
hold off

En la gráfica anterior, podemos apreciar como en la parte inferior de la placa el más caliente, y cuanto más nos acercamos a la líneas de temperatura circulares la cual vemos en el centro de la parte superior, más se enfría esta. Esto podría coincidir con una placa de metal calentado al cual le está cayendo un chorro de agua (o viento frío) en el centro de dicho círculo.

3 Gradiente de la temperatura.

Representación del gradiente de la temperatura.

El código adjunto describe la forma en que hemos calculado el gradiente del campo escalar temperatura. A la derecha vemos representado dicho campo en sus curvas de nivel con el gradiente siendo perfectamente ortogonal a las curvas isotermas.

%---Divergencia del campo de temperaturas---
r=linspace(0,10,50);
t=linspace(0,2,10);
[x,y]= meshgrid(r,t);      %hasta aquí es lo mismo que el resto de apartados.

u=2.*(x-6) ;      %componente i
v=20.*(y-3/2);   %componente j

T=((x-6).^2)+(10.*(y-3/2)).^2;

hold on       %Campo escalar de temperatura            
 contour(x,y,T,30)             %Líneas de nivel
colorbar                      %Mostramos las escala
axis([-0.5,10.5,-0.5,2.5])
w=quiver(x,y,u,v);
axis([-0.5,10.5,-0.5,2.5])
title ('Gradiente de la Temperatura','Fontsize',18);
set(w,'maxheadsize',0.33)
xlabel('Eje X')
ylabel('Eje Y')
hold off

Representación del gradiente de la temperatura.

Si aumentamos la imagen se puede apreciar que las flechas que representan el gradiente son ortogonales a sus respectivas líneas de campo. Esto ha de cumplirse por definición del gradiente de un campo y muestra la dirección de máximo aumento de temperatura en cada punto.

4 Campo de deformaciones para el instante inicial.

Representación del mallado.

%----Campo de deformaciones____
r=linspace(0,10,50);
t=linspace(0,2,10);
%mallado
[x,y]= meshgrid(r,t);      
          
Z=zeros(10,50);
hold on
surf(x,y,Z)
 u=1.*x ;      %componente i
 v=(1/3).*sin(pi/3.*y) ;   %componente j
w=quiver(x,y,u,v,'r');
axis([-0.5,10.5,-0.5,2.5])
title ('Gradiente de la Temperatura','Fontsize',18);
set(w,'maxheadsize',0.33)
xlabel('Eje X')
ylabel('Eje Y')
view(2)
hold off

Este gráfico muestra la dirección e intensidad de las deformaciones, en el podremos apreciar zonas de convergencia y divergencia.

5 Placa antes y después del desplazamiento

Al analizar el campo de deformaciones del apartado anterior, podemos suponer que se formaran ondas a lo largo de la placa. Es por eso que para ver bien la variación tendremos que comprobar en una representación tridimensional.

Placa antes y después del desplazamiento

%Posición de la placa deformada:
  rx=Mx+(2/5).*sin(Mx);
  ry=0.*My; 
  %Representación de la placa antes del desplazamiento: 
  subplot(1,2,1)
  surf(Mx,My,0*Mx)
  title('Posición original')
  axis equal
  axis([-3,3,-1,5])
 view(2) 
 xlabel('Eje X')
 ylabel('Eje Y') 
 zlabel('Eje Z')
 %Representación de la placa después del desplazamiento:
 subplot(1,2,2)
 surf(rx,ry,0*rx)
 title('Después del desplazamiento')
 axis equal
 axis([-3,3,-1,5])
 xlabel('Eje X')
 ylabel('Eje Y')  
 zlabel('Eje Z') 
 view(2)

6 Divergencia del campo vectorial sobre la capa.

Para el cálculo de la divergencia en coordenadas cartesianas, se procede a derivar el campo de desplazamientos hallado anteriormente:

[math]\ \vec u (x, y, z)= \frac{2}{5} \sin (x) \vec i [/math]

y por tanto:

[math]\ \nabla \cdot \vec u =\frac{\partial \vec u_1}{\partial x} + \frac{\partial \vec u_2}{\partial y} + \frac{\partial \vec u_3}{\partial z} = \frac{\partial (\frac{2}{5}\sin (x))}{\partial x} + \frac{\partial 0}{\partial y} + \frac{\partial 0}{\partial z} = \frac{2}{5} \cos(x) [/math]

Representación de la divergencia del campo vectorial sobre la capa.

clear
clc
h=0.2;
x=0:h:10;           
y=0:h:2;
%Creamos el mallado
[Mx,My]=meshgrid(x,y);  
%la z es cero por ser una placa sin espesor
z=zeros(size(Mx)); 
mesh(Mx,My,z)         
%Pintamos las curvas de nivel correspondientes a la temperatura
g=(2/5)*cos(Mx);        
pcolor(Mx,My,g)
%para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
shading flat             
contour(Mx,My,t,30,'r') 
colorbar
hold on
%dibujamos                 
quiver3(Mx,My,z,u,v,w)             
axis([-0.5,10.5,-0.5,2.5])
view(2)
hold off

7 Calcular |∇ × [math]\vec{u}[/math] | en todos los puntos del solido en t = 0

Sea [math]|∇ × \vec{u}|[/math] el rotacional de un campo de desplazamientos [math] \vec u = (ux, uy, uz)[/math] expresado en un sistema de coordenadas de vectores unitarios ortogonales [math]\vec{a}, \vec{b}, \vec{c}[/math], aplicado a una superficie bidimensional expresado en dicho sistema de coordenadas, se cumple que:

[math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{u} & \vec{v} &\vec{w} \\ d/da & d/db & d/dc\\ ux & uy & uz \end{vmatrix}[/math]

Por ello, para calcular el rotacional de un campo de desplazamientos, se aplica la fórmula del producto vectorial en los ejes del sistema de coordenadas.

Para el caso del sistema de coordenadas cartesiano, con ejes [math][ x, y, z ][/math], y vectores respectivos [math]\vec{i}, \vec{j}, \vec{k}[/math], se procede a calcular el rotacional.

Usando el campo [math]\vec{u} = (\vec{ux}, \vec{uy}, \vec{uz}) = (2/5sin(x), 0, 0)[/math] previo, procedemos a hacer los cálculos:

[math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{i} & \vec{j} &\vec{k} \\ d/dx & d/dy & d/dz\\ ux & uy & uz\end{vmatrix}= \begin{vmatrix}\vec{i} & \vec{j} & \vec{k} \\ d/dx & d/dy & d/dz \\ 2/5sin(x) & 0 & 0\end{vmatrix}=(0-0)\vec{i}+(0-0)\vec{j}+(0-0)\vec{k}=0 [/math]

Por lo tanto, se trata de un campo conservativo es decir, el rotacional en todos los puntos de la placa es nulo, lo que implica:

[math]\ \left | \bigtriangledown \times \vec{u} \right | = 0 [/math]

Dirección de los vectores de deformacion.

Además este resultado tiene un significado gráfico. Al representar los vectores de deformación apreciamos que no aparece ninguna rotación, en su lugar solo aparecen zonas de divergencia y convergencia (divergencia negativa).

8 Tensor de deformaciones

Definiendo [math]ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}[/math], la parte simétrica del tensor gradiente de [math]\vec u[/math] conocido como tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones a través de la fórmula [math]σ=λ\nabla·\vec u 1+2με[/math].

[math]ε(\vec u)=\nabla \vec u= \begin{pmatrix} \frac{2cosx}{5} & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]σ=λ\nabla·\vec u 1+2με = λ\frac{2cosx}{5}·\begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} + 2μ\begin{pmatrix}\frac{2cosx}{5} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = 2λ\frac{cosx}{5}·I+2μ·\frac{2cosx}{5}\begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} = [/math]

[math]= \frac{2cosx}{5} (\begin{pmatrix} λ & 0 & 0\\ 0 & λ & 0 \\ 0 & 0 & λ \end{pmatrix} + \begin{pmatrix} 2μ & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}) = \frac{2cosx}{5}\begin{pmatrix} 2μ+λ & 0 & 0\\ 0 & λ & 0 \\ 0 & 0 & λ \end{pmatrix}→{λ=μ=1}→\frac{2cosx}{5}\begin{pmatrix} 3 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math]

Tomando [math]λ = µ = 1[/math], calculamos las tensiones normales en la dirección que marca el eje [math]\vec i[/math], es decir [math]\vec i · σ ·\vec i[/math], las tensiones normales en la dirección que marca el eje
[math]\vec j[/math], es decir [math]\vec j · σ · \vec j [/math] y las correspondientes al eje [math]\vec k[/math], es decir [math]\vec k · σ ·\vec k[/math].

Proseguimos con el cálculo de las tensiones normales:

• [math]\vec i · σ ·\vec i = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}·\begin{pmatrix} \frac{6cosx}{5} & 0 & 0\\ 0 & \frac{2cosx}{5} & 0 \\ 0 & 0 & \frac{2cosx}{5} \end{pmatrix}·\begin{pmatrix} 1 \\0 \\ 0 \end{pmatrix} = \frac{6cosx}{5}[/math]

• [math]\vec j · σ ·\vec j = \begin{pmatrix} 0 & 1 & 0 \end{pmatrix}·\begin{pmatrix} \frac{6cosx}{5} & 0 & 0\\ 0 & \frac{2cosx}{5} & 0 \\ 0 & 0 & \frac{2cosx}{5} \end{pmatrix}·\begin{pmatrix} 0 \\1 \\ 0 \end{pmatrix} = \frac{2cosx}{5}[/math]

• [math]\vec k · σ ·\vec k = \begin{pmatrix} 0 & 0 & 1 \end{pmatrix}·\begin{pmatrix} \frac{6cosx}{5} & 0 & 0\\ 0 & \frac{2cosx}{5} & 0 \\ 0 & 0 & \frac{2cosx}{5} \end{pmatrix}·\begin{pmatrix} 0 \\0 \\ 1 \end{pmatrix} = \frac{2cosx}{5}[/math]

Tensor de deformaciones.

%____Tensor de deformaciones____
x=linspace(0,10,50);      
y=linspace(0,2,10);     
[Mx,My]= meshgrid(x,y); 
Ti=(6/5).*cos(Mx);       % Tension normal en i
Tj=(2/5).*cos(Mx);       % Tension normal en j
Tk=(2/5).*cos(Mx);       % Tension normal en k
figure
subplot(1,3,1)           
pcolor(Mx,My,Ti)
colorbar
shading flat
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
axis([-0.5,10.5,-0.5,2.5]);
title('Tensiones normales en i')  
subplot(1,3,2) 
pcolor(Mx,My,Tj)
colorbar
shading flat
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
axis([-0.5,10.5,-0.5,2.5]);
title('Tensiones normales en j')
subplot(1,3,3)
pcolor(Mx,My,Tk)
colorbar
shading flat
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
axis([-0.5,10.5,-0.5,2.5]);
title('Tensiones normales en k')

9 Tensiones tangenciales

Sea [math]\vec{T}[/math] la tensión tangencial respecto al plano ortogonal a [math]\vec{i}[/math], que viene dado por la siguiente expresión:

[math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math]

Particularizando para nuestros datos obtenemos:

[math]\frac{6}{5}\cos (x)\vec{i}-\frac{6}{5}\cos (x)\vec{i}=\vec{0}[/math]

Por lo tanto la tensión tangencial respecto al plano que contiene a los vectores [math]\vec{j}[/math] y [math]\vec{k}[/math] es nula.

10 Tensión de Von Mises

La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión, la teoría expone que un material dúctil comienza a ceder en una ubicación cuando la tensión de Von Mises es igual al límite de tensión. En la mayoría de los casos, el límite elástico se utiliza como el límite de tensión. En este caso podemos observar que el punto máximo de tensión es 0,9791.
La tensión de Von Mises puede calcularse fácilmente a partir de las tensiones principales del tensor tensión en un punto de un sólido deformable, mediante la expresión:

[math]σ=\sqrt{\frac {(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}{2}}[/math]

Donde [math]σ_1,σ_2,σ_3[/math] son los autovalores de [math]σ=\frac{2cosx}{5}\begin{pmatrix} 3 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math].

Tensión de Von Mises en 3D.

x=linspace(0,10,50);                     % Vector x con valores entre 0 y 10    
y=linspace(0,2,50);                      % Vector x con valores entre 0 y 2      
[X,Y]=meshgrid(x,y);                     % Obtención del mallado      
% Definimos nuestra función                                               
fx=inline('(6/5)*cos(x)','x','y');        
fy=inline('(2/5)*cos(x)','x','y');
fz=inline('(2/5)*cos(x)','x','y');
% Inicio del bucle
for i= 1:length(x)                              
   for j= 1:length(y)                             
% 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
% Creación de un vector con las componentes [xx, yy, zz]
v=[xx yy zz];                            
M=diag(v);                               % Diagonalización del vector       
autov=eig(M);                            % Obtención de los autovalores      

% Cálculo de la tensión de Von Mises para cada componente                                                
VonM=sqrt(((autov(1)-autov(2))^2+((autov(2)-autov(3))^2+((autov(3)-autov(1))^2))*1/2));  
Z(i,j)=VonM;
   end
end
% Fin del bucle                                             

subplot(2,1,1)
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

subplot(2,1,2)
plot(X,Z,'*r')                                  % Representación superficial con líneas entrecruzadas
title('TENSIÓN DE VON MISES EN 2D')             % Título del gráfico
xlabel('Eje X')                                 % Título del eje x
ylabel('Eje Z')                                 % Título del eje y

maximo=max(max(Z))                           % Obtención del valor máximo

11 Velocidad de propagación v en términos de Lamé

Queremos calcular la velocidad de propagación de las ondas v en términos de las constantes de Lamé , por lo cual tendremos que suponer [math]\vec{F}[/math]=0.

El desplazamiento en la placa viene dado por las fuerzas [math]\vec{F}[/math] que actúan sobre dicha placa, las cuales se pueden aproximar con la siguiente función de elasticidad:

[math]\vec{F} = \frac{d^2\vec{u}}{dt^2} - \bigtriangledown \cdot \sigma[/math]

definimos en nuestro caso que:

[math]\vec{u}[/math] = [math]\frac{2}{5}\sin(x-vt)\vec{i}[/math]

siendo v la velocidad de propagación de las deformaciones y como hemos calculado en el apartado 8, [math]\sigma[/math] será :

[math]\sigma = \begin{pmatrix} \frac{6}{5}cos(x-vt) & 0 & 0 \\ 0 &\frac{2}{5}cos(x-vt) & 0 \\ 0 & 0 & \frac{2}{5}cos(x-vt)\end{pmatrix}[/math]

De esta forma la divergencia [math]\nabla \cdot \sigma[/math] nos quedaría de la siguiente forma:

[math] \nabla \cdot \sigma = \nabla \cdot (\frac{6}{5}cos(x-vt)\vec{i}) + \nabla \cdot (\frac{2}{5}cos(x-vt)\vec{j}) + \nabla \cdot (\frac{2}{5}cos(x-vt)\vec{k}) = -\frac{6}{5}sin(x-vt)-\frac{2}{5}sin(x-vt)-\frac{2}{5}sin(x-vt)=-2 sin(x-vt)[/math]

Así nos quedaría que:

[math]\vec{F} = \frac{d^2\vec{u}}{dt^2} -2 sin(x-vt) = 0[/math] [math]\vec{F} = \frac{d}{dt}(-\frac{2}{5}v cos(x-vt)) -2 sin(x-vt) = -\frac{2}{5}v^2 sin(x-vt) -2 sin(x-vt) = 0[/math]

Para que se cumpla que para todo [math]x[/math] y todo [math]t[/math], [math]F = 0[/math] solo es posible de dos maneras, siendo v un número imaginario que cumpla que [math]-\frac{2}{5}v^{2}-2=0[/math] o un valor que dependa de [math]x[/math] y [math]t[/math], el cual cumpla que [math]sin(x-vt) = 0[/math];

• [math]sin(x-vt) = 0[/math] la solución de esta ecuación es [math]v=\frac{x}{t}[/math] resultado imposible para [math]t=0[/math].

• [math]-\frac{2}{5}v^{2}-2=0[/math] esta ecuación nos presenta la siguiente solución imaginaria: [math]v=\sqrt{5} i[/math]

12 Módulo de desplazamiento en dirección [math]\vec{j}[/math]

La onda es longitudinal, es decir, que se desplaza en el mismo sentido que la amplitud, por tanto no habrá desplazamiento en otro eje distinto al eje en el que se transmite dicha onda, de aquí la demostración:

Adoptando el campo de desplazamientos: [math]\vec u (x,y,t)=\frac{2}{5}\vec i sen(x-vt)[/math]

Fijando el punto [math](x,y)=(1/2,1)[/math] se tiene: [math]\vec u (1/2,1,t)=\frac{2}{5}\vec i sen(1/2-vt)[/math] donde el parámetro [math]t[/math], tiempo, pertence al intervalo [math][0,10][/math]

En la dirección de [math]\vec j [/math] se obtiene que: [math]\vec u \cdot \vec j = \frac{2}{5}sen(1/2-vt)\vec i \cdot \vec j = 0[/math] como se esperaba por tratarse de una onda longitudinal.

En la dirección de [math]\vec i [/math] se obtiene que: [math]\vec u \cdot \vec i = \frac{2}{5}sen(1/2-vt)\vec i \cdot \vec i = \frac{2}{5}sen(1/2-vt)[/math]

Es decir, al ser una onda longitudinal en la que la dirección de propagación es la misma que la amplitud, se producen desplazamientos en el eje [math]\vec{i}[/math] y no en la dirección [math]\vec{j}[/math], demostrando así lo dicho arriba.