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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad.
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Erick Morales Pruna Hugo Sacristán de Agustín Jaime Villalba Guerrero Ángel Matín Cruz
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Visualización de campos escalares y vectoriales en elasticidad. Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región [math](x, y) ∈ [-1, 1]×[0, 12][/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) = log(1+x^2)+log(1+(y-4)^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)= x \vec{i} + y \vec{j} [/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_{d}}(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 ondulatorio de los puntos de la misma dado por el vector:

[math]\vec{u}(x, y, t)=\vec{a}sin(k\pi(\vec{d}·\vec{r_{0}}(x,y)-vt)),[/math]donde [math]\vec{a}[/math] se conoce como amplitud, k>0 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.

La variable t representa el tiempo que congelaremos en t=0 en los primeros apartados de este trabajo de manera que supondremos, para los primeros apartados,

[math]\vec{u}(x, y, t)=\vec{a}sin(k\pi(\vec{d}·\vec{r_{0}}(x,y))).[/math]

Supondremos que se trata de una onda transversal en la que la dirección de propagación es ortogonal a la amplitud. Tomaremos en particular:

[math]\vec{a}= 1/3 \vec{i}, k=1, \vec{d}= 1/3 \vec{j}[/math]

1 Representación de la placa rectangular plana.

Dibujar un mallado que represente los puntos interiores del sólido. Tomar los ejes (comando axis) en el rectángulo [math](x, y) ∈ [−1; 1] × [0;12][/math] y como paso de muestreo [math]h = 2/10[/math] para las variables [math]x[/math] e [math]y[/math].

Figura 1: Mallado de la placa.

% Paso de muestreo h para las variables x e y.
h=2/10;
x=[-1:h:1];
y=[0:h:12];
% Mallado con las matrices Mx e My.
[Mx,My]=meshgrid(x,y);
mesh(Mx,My,0*My);
% Ejes
axis([-5,5,-0.5,12.5]);
% Escribimos el titulo del gráfico y los nombres de los ejes.
title('Mallado de la placa');
xlabel('Eje X');
ylabel('Eje Y');
% Con el comando view(2), visualizamos el mallado en 2 dimensiones.
view(2);

2 Representación de las curvas de temperatura.

Dibujar las curvas de nivel de la temperatura (comando contour) y decidir en qué punto la temperatura es máxima a partir de la gráfica.

Primeramente calcularemos el gradiente de la temperatura con la siguiente formula:

[math]∇T=\frac{\partial T}{\partial x}\vec{i} + \frac{\partial T}{\partial y}\vec{j} = \frac{2x}{1+x^2}\vec{i}+\frac{2y-8}{1+(y-4)^2}\vec{j}[/math].

Tal y como podemos observar en la dirección de las flechas de la figura 2 estas son siempre perpendiculares a las curvas de nivel, esto se debe a que el gradiente de un vector nos indica la dirección de máximo crecimiento en cada punto del vector.

Para hallar cual es la máxima temperatura utilizaremos el comando [math]max(max(T))[/math], dicha temperatura se alcanzará en dos puntos, en el (x=1,y=12) y en el (x=-1,y=12) tal y como podemos obsrvar en la gráfica.

Figura 2: Representación curvas de nivel de la temperatura.

h = 2/10;
x = -1:h:1;
y = 0:h:12;
% Creación del mallado
[Mx,My]= meshgrid(x,y);
% Función temperatura
T =log(1+(Mx.^2))+log(1+(My-4).^2);
% Se define el rango de visión de la gráfica.
axis([-1,1,0,12]);
% Representación de la temperatura y las curvas de nivel
subplot(1,2,1);
mesh(Mx,My,T);
subplot(1,2,2);
contour(Mx,My,T,20);
colorbar
hold on
x=-1:h:1;     
y=0:h:12;           
[Mx,My]=meshgrid(x,y); 
figure(1)
% Gradiente de T
fx=(2.*Mx)./(1+(Mx.^2));               
fy=((2.*My)-8)./(1+(My-4).^2);                
quiver(Mx,My,fx,fy)     
axis([-1,1,0,12])     
view(2)

3 Ley de Fourier

De acuerdo a la Ley de Fourier la energía calorífica [math]\vec{Q}[/math] viaja de acuerdo a la fórmula : [math]\vec{Q}=-k*∇T,[/math] donde k es la constante de conductividad térmice de la placa que supondremos k=1. Calcular [math]\vec{Q}[/math] y dibujarlo como campo vectorial.

Finalmente obtenemos que [math]\vec{Q}[/math] es igual a:

[math] \vec{Q} = -\frac{2x}{1+x^2}\vec{i}-\frac{2y+8}{1+(y-4)^2}\vec{j}[/math].

Figura 3: Energía calorífica.

h = 2/10;
x =-1:h:1;
y = 0:h:12;
% Creación del mallado
[Mx,My]= meshgrid(x,y);
% Función temperatura
T =log(1+(Mx.^2))+log(1+(My-4).^2);
% Se define el rango de visión de la gráfica.
axis([-1,1,5,12]);
% Gradiente de T
fx=(2.*Mx)./(1+(Mx.^2));               
fy=((2.*My)-8)./(1+(My-4).^2);  
qx=-1.*fx;
qy=-1.*fy;
% Título
title('Energía calorífica');
% Representación de la temperatura y las curvas de nivel
hold on
quiver(Mx,My,qx,qy)
contour(Mx,My,T,20);
colorbar

4 Representación del campo de vectores en t=0.

Dibujar el campo de vectores en los puntos del mallado del sólido, en t = 0.

Figura 4: Campo de desplazamientos.

h = 2/10;
x =-1:h:1;
y = 0:h:12;
%Creación del mallado
[Mx,My]= meshgrid(x,y);
%Componentes en la dirección de i y de j del campo de desplazamiento
ux=(1/3).*sin((pi/3).*My) ;
uy=0.*My ;
figure
%dibujo del mallado 
mesh(Mx,My,0*Mx)
hold on 
%campo de desplazamientos
quiver(Mx,My,ux,uy,'k')
axis([-2,2,-2,15])
view(2)
title('Campo de desplazamientos')
hold off

5 Representación del desplazamiento del sólido.

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec{u}[/math] (en t = 0). Dibujar ambos en la misma figura usando el comando subplot.

Figura 5: Sólido antes y despues del desplazamiento.

h= 2/10;
x=-1:h:1;
y= 0:h:12;
%Creación de matriz x e y
[Mx,My]=meshgrid(x,y);
%posicion final 
rx=((1/3).*sin((pi/3).*My))+Mx;
ry=(0.*My)+My; 
%representacion de la superficie antes del desplazamiento  
subplot(1,2,1)
surf(Mx,My,0*Mx)
title('Antes del desplazamiento')
axis([-1.5,1.5,-0,13])
view(2) 
xlabel('x')
ylabel('y') 
zlabel('z')
%representacion de la superficie después del desplazamiento
subplot(1,2,2) 
surf(rx,ry,0*rx)
title('Después del desplazamiento')
axis([-1.5,1.5,0,13])
view(2);
xlabel('x')
ylabel('y') 
zlabel('z')

6 Estudio analítico de la divergencia

Dibujar [math]∇·\vec{u}[/math] en [math]t=0[/math]. Determinar analíticamente los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. La divergencia es una medida del cambio de volumen local debido al desplazamiento. ¿Se puede apreciar esto en la gráfica?

Primeramente realizamos la divergencia aplicando la siguiente fórmula:

[math]∇·\vec{u} = \frac{∂}{∂x}(u_1)+\frac{∂}{∂y}(u_2)=0.[/math].

Sin embargo, como el vector [math]\vec{u}[/math] es [math]\vec{u}=\frac{1}{3}·sen(\frac{\pi y}{3})\vec{i}[/math], la divergencia es 0.

La razón por la cual la divergencia es 0 se debe a que el campo [math]\vec{u}[/math] es senoidal.

7 Cálculo y representación del rotacional de u.

Calcular [math]|∇ × \vec{u}|[/math] en todos los puntos del sólido en [math]t = 0[/math] y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

[math]∇×\vec u(x,y,z) = \begin{vmatrix} \vec{e_i} & \vec{e_j} & \vec{e_k} \\ \frac{∂}{∂x} & \frac{∂}{∂y} & \frac{∂}{∂z} \\ \frac{1}{3}·sen(\frac{\pi y}{3}) & 0 & 0\end{vmatrix} = -\frac{\pi}{9}cos(\frac{\pi y}{3})\vec{e_k}[/math]

Una vez hallado el rotacional, hallamos el módulo: [math]\frac{\pi}{9}cos(\frac{\pi y}{3})[/math]

Tal y como podemos observar en la gráfica el valor del rotacional oscila a lo largo de la gráfica, alcanzo el valor más alto del rotacional en dos puntos:

- [math]P_1(x,y,z) = (-1\leq[/math]x[math]\leq1,y=3,z=0.35).[/math]

- [math]P_2(x,y,z) = (-1\leq[/math]x[math]\leq1,y=9,z=0.35).[/math]

Figura 6: Módulo del rotacional.

h= 2/10;
x=-1:h:1;
y= 0:h:12;
%Creación de matriz x e y
[Mx,My]=meshgrid(x,y); 
%Módulo del rotacional
rot = (-pi/9).*cos((pi/3).*My);
%Representación gráfica del rotacional
surf(Mx,My,rot)
shading flat
axis equal
colorbar
view(3);
axis([-1.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('X');
ylabel('Y');
zlabel('z');

8 Tensor de deformaciones.

Definamos [math]ϵ(\vec{u}) = (∇\vec{u} + ∇\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 [math]σ_{ij}[/math] a través de la fórmula: [math]σ = λ∇ · \vec{u} 1 + 2µϵ[/math]

donde [math]1[/math] es el tensor identidad en el conjunto de vectores libres del espacio [math]R^3[/math] y [math]λ[/math], [math]µ[/math] 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 [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. Tomando [math]λ = µ = 1[/math], dibujar 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] (dibujar las que no son nulas).

Primeramente calcularemos [math]ϵ(\vec{u})[/math] = [math]Ԑ(\vec{u}) = \begin{pmatrix} 0 & \frac{\pi}{18}cos(\frac{\pi y}{3}) & 0 \\ \frac{\pi}{18}cos(\frac{\pi y}{3}) & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

Después sabiendo que [math]∇ · \vec{u}=0[/math], λ=1 y µ=1, aplicamos la siguiente fórmula para hallar el tensor de tensiones, obteniendo lo siguiente:

[math]σ = λ∇ · \vec{u} 1 + 2µϵ = 2·Ԑ(\vec{u})= \begin{pmatrix} 0 & \frac{\pi}{9}cos(\frac{\pi y}{3}) & 0 \\ \frac{\pi}{9}cos(\frac{\pi y}{3}) & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

Una vez hallado el tensor de tensores hallaremos las tensiones normales en las direcciones que marca el eje [math]\vec{i},\vec{j},\vec{k}[/math], obteniendo:

[math]\vec{i}· σ · \vec{i} = \vec{j}· σ · \vec{j} = \vec{k}· σ · \vec{k} = 0 [/math].

Debido a que en todas las direcciones tienen un valor nulo, no es posible su representación.

9 Tensiones tangenciales.

Calcular las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math], es decir [math]|σ ·\vec{i} − (\vec{i} · σ ·\vec{i})\vec{i}|[/math], en [math]t = 0[/math]. Dibujar sólo las que no son nulas.

[math]|σ·\vec{i} − (\vec{i}·σ·\vec{i})\vec{i}|= |σ·\vec{i}| =|\begin{pmatrix} 0 & \frac{\pi}{9}cos(\frac{\pi y}{3}) & 0 \\ \frac{\pi}{9}cos(\frac{\pi y}{3}) & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\begin{pmatrix} 1\\0\\0 \end{pmatrix}|= |\begin{pmatrix} 0\\ \frac{\pi}{9}cos(\frac{\pi y}{3}) \\0 \end{pmatrix}|[/math]

Finalmente la tensión tangencial respecto al plano ortogonal a [math]\vec{i}[/math] es [math]\frac{\pi}{9}cos(\frac{\pi y}{3})[/math].

En la siguiente figura podemos obstervar la tensión tangencial respecto al plano ortogonal a [math] \vec{i}[/math].

Figura 7: Tensiones tangenciales.

h= 2/10;
x=-1:h:1;
y= 0:h:12;
%Creación de matriz x e y
[Mx,My]=meshgrid(x,y); 
%Tensión tangencial en cada punto
tn=(pi/9).*cos((pi/3).*My);
%Representación gráfica
quiver(Mx,My,tn,tn.*0);
axis([-1.5,1.5,-1.5,13.5]);
title('Tension tangencial');
xlabel('X');
ylabel('Y');

10 Tensión de Von Mises

La tensión de Von Mises se define por la fórmula [math]σ_{VM}=\sqrt{\frac{(σ_{1}-σ_{2})^2+(σ_{2}-σ_{3})^2+(σ_{3}-σ_{1})^2}{2}}[/math]

donde [math]σ_{1}[/math], [math]σ_{2}[/math] y [math]σ_{3}[/math] son los autovalores de [math]σ[/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 (y no elástico puro). Pintar la tensión de Von Mises y señalar en qué punto se alcanza el mayor valor. (Para calcular autovalores con OCTAVE/MatLab usar el comando eig.m)

Tal y como podemos obsvervar en la gráfica y con la ayuda de MATLAB encontramos los puntos en los que se alcanza la máximo tensión de Von Mises, siendo estos puntos los siguientes, [math]y=0, y=3, y=6, y=9[/math] e [math] y=12. [/math]

Figura 8: Tensión de Von Mises.

h= 2/10;
x=-1:h:1;
y= 0:h:12;
%Creación de matriz x e y
[Mx,My]=meshgrid(x,y); 
%definimos la función de Von mises. t1,t2,t3 son las tensiones principales
VonMises=inline('(((t1-t2)^2+(t2-t3)^2+(t3-t1)^2)/2)^(1/2)','t1','t2','t3');
[f,c]=size(Mx);
%asignamos a la matriz MVonM los valores de la tensión de Von Mises en cada punto
  for i=1:f
   for j=1:c
    deformaciones=[[0;(pi/9).*cos((pi/3).*My(i,j));0],[(pi/9).*cos((pi/3).*My(i,j));0;0],[0;0;0]];
    sigmas=eig(deformaciones);
    t1=sigmas(1,1);
   t2=sigmas(2,1);
   t3=sigmas(3,1);
   Mvon(i,j)=VonMises(t1,t2,t3);
  end
  end
%Representación gráfica
surf(Mx,My,Mvon)
axis([-1.5,1.5,-0.5,12.5]);
shading interp
axis equal
title('Tensión de Von Mises');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
view(3);
colorbar

11 Campo de fuerzas que actúa sobre la placa

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

[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ[/math]

donde [math]∇ · σ[/math] es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz [math]σ[/math].

Calcular la velocidad de propagación de las ondas [math]v[/math] en términos de las constantes de Lamé, suponiendo que [math]\vec{F} = 0[/math]. Si la onda fuera longitudinal, es decir, tomando [math]\vec{a} = 1/3\vec{j}[/math], ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

Para este apartado tendremos que utilizar el siguiente vector [math] \vec{u}[/math]:

[math]\vec{u}(x,y,t)=\vec{a}·sin(\pi k(d·\vec{r0}(x, y)−vt))[/math]

siendo [math]\vec{a}=\frac{1}{3}\vec{i}, k=1, \vec{d}=\frac{1}{3}\vec{j}.[/math]

Sustituyendo obtenemos que [math]\vec{u}(x,y,t)=\frac{1}{3}·sin(\frac{\pi y}{3}-\pi vt)\vec{i})[/math]

Para hallar el campo de fuerzas [math]\vec{F}[/math], deberemos redefinir los siguientes parametros:

[math]∇·\vec{u} = \frac{∂}{∂x}(u_1)+\frac{∂}{∂y}(u_2)=0.[/math]

[math]ϵ(\vec{u})[/math] = [math]Ԑ(\vec{u}) = \begin{pmatrix} 0 & \frac{\pi}{18}cos(\frac{\pi y}{3}-\pi vt)) & 0 \\ \frac{\pi}{18}cos(\frac{\pi y}{3}-\pi vt) & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]σ = λ∇ · \vec{u} 1 + 2µϵ = 2·Ԑ(\vec{u})= \begin{pmatrix} 0 & \frac{\pi}{9}cos(\frac{\pi y}{3} - \pi vt) & 0 \\ \frac{\pi}{9}cos(\frac{\pi y}{3} - \pi vt) & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]∇ · σ[/math] queda [math]∇ · σ = \frac{\pi ^2}{27} · sen(\frac{\pi y}{3}-\pi vt)\vec{i}[/math]

[math]\frac{∂\vec{u}}{∂t}= -\frac{\pi v}{3}·cos(\frac{\pi y}{3}-\pi vt)\vec{i} [/math]

[math]\frac{∂^2\vec{u}}{∂t^2}= -\frac{\pi ^2 v^2}{3}·sin(\frac{\pi y}{3}-\pi vt)\vec{i}[/math]

Una vez hecho esto podemos aplicar la ecuación de la elasticidad lineal: [math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ[/math], con la cual metiendo los datos obtenidos, suponiendo que [math]\vec{F}=0[/math] y despejando la velocidad obtendremos una velocidad de propagación de las ondas de, [math]\vec{v}=\frac{1}{3}\vec{i}[/math].

Si la onda longitudinal fuera [math]\vec{a}=\frac{1}{3}\vec{j}[/math], repetiríamos el mismo proceso, obteniendo una velocidad de onda longitudinal de [math]\vec{v}=\frac{1}{\sqrt{3}}\vec{j}.[/math]

Tal y como podemos observar en los cálculos, estando en el mismo medio las ondas transversales y longitudinales tienen una velocidad diferente.

12 Módulo de desplazamiento transversal

Fijado ahora el punto [math](1/2, 1)[/math], calcular el módulo del desplazamiento transversal (dirección [math]\vec{i}[/math]) a lo largo del tiempo en el intervalo [math]t ∈ [0, 10][/math]. Dibujar la función que a cada t le asocia dicho desplazamiento.

Sustituyendo el punto P en el vector [math]\vec{u}[/math] obtenemos lo siguiente: [math]\vec{u}(P)=\frac{1}{3}·sen(\frac{\pi y}{3}-\pi vt)\vec{i}[/math].

Para hallar el módulo del desplazamiento transversal hacemos estas operaciones: [math]\vec{u}·\vec{i}=\frac{1}{3}·sen(\frac{\pi}{3}-\frac{\pi}{3}t)[/math], siendo este un desplazamiento horizontal. [math]\vec{u}·\vec{j}=0[/math], tal y como podíamos esperar por tratarse de una onda transversal, no habiendo desplazamiento vertical.

Figura 9: Representación de la función desplazamiento.

h= 2/10;
t=0:h:10;
f=@(t) (1/3)*sin((pi/3)-((pi/3).*t));
desplazamiento=f(t);
plot(t,desplazamiento)
title ('Representación de la función desplazamiento')