Trabajo 2. Visualización de Campos Escalares. Grupo 29-C
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales en elasticidad. Grupo 29-C |
| Asignatura | Teoría de Campos |
| Curso | 2022-23 |
| Autores | Carlos Eduardo Sayalero Gómez. Caroline Vidal Rainho. Diego Andrés Vásquez Peñaloza. Sergio Alejandro Pareja barbarito. |
| 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 [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:
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
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 [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.
En este caso, tendremos los siguientes valores: [math] \vec{a}=\frac{2}{5}\vec{i}[/math], [math]\vec{d}=\vec{i}[/math] y [math]k=1 [/math] , con los cuales obtenemos el siguiente vector [math]\vec{u}(x,y,t)[/math]
1. Mallado de la placa.
Dibujar un mallado que represente los puntos interiores del sólido. 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 = 2/10[/math] para las variables [math]X[/math] e [math]Y[/math].
% Dibujar un mallado que represente los puntos del sólido
% Definimos la malla a partir de los vectores en coord. cartesianas:
x=0:0.2:10;
y=0:0.2:2;
% Realizamos el mallado:
[xx,yy]=meshgrid(x,y);
% Dibujamos la placa:
mesh(xx,yy,0*xx)
% Ajustamos los ejes para que tengan la misma escala:
axis equal
% Ajustamos los ejes según el enunciado:
axis([-0.5 10.5 -0.5 2.5])
% Para ver la placa desde arriba:
view(2)
2. Temperatura de la placa.
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.
En este apartado se nos pide representar las curvas de nivel de la temperatura de la placa y, a partir de los gráficos, decir en qué punto se encuentra la temperatura máxima.
% Dibujar las curvas de nivel
% Discretizamos variables
x=0:0.2:10;
y=0:0.2:2;
% Realizamos el mallado:
[xx,yy]=meshgrid(x,y);
% Calculamos la temperatura para cada punto del mallado
Temp=(xx-6).^2+(10*(yy-3/2)).^2;
% Dibujamos las curvas de nivel en la placa
hold on
mesh(xx,yy,0*xx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
% Para ver la placa desde arriba:
view(2)
[M,c]=contour(xx,yy,Temp,[0 1.1 2 5 10 20 30 40 50 100 150 200],'ShowText','on')
hold off
Como se puede apreciar en la figura, en la parte inferior del sólido, las temperaturas alcanzadas son mayores que en la parte superior. Esto se puede dar, entre otras cosas, porque el sólido está siendo calentado desde su parte inferior.
3. Gradiente de la temperatura.
Calcular [math]∇T[/math] y dibujarlo como campo vectorial. Observar gráficamente que [math]∇T[/math] es ortogonal a dichas curvas.
Realizaremos el calculo del gradiente:
% Representación del gradiente de temperatura
x=0:0.2:10;
y=0:0.2:2;
% Realizamos el mallado:
[xx,yy]=meshgrid(x,y);
% Calculamos la temperatura para cada punto del mallado
Temp=(xx-6).^2+(10*(yy-3/2)).^2;
% Calculamos las derivadas parciales de T
GradTx=2.*(xx-6);
GradTy=200.*(yy-3/2);
% Dibujamos el gradiente
mesh (xx,yy,0*xx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
contour(xx,yy,Temp,[0 1.1 5 20 50 150 200 240 250 260]);
quiver(xx,yy,GradTx,GradTy,1,'r');
En esta ocasión, es necesario incorporar otra imagen, la cual hace zoom en la parte superior derecha de la imagen original. Esto se hace con la intención de apreciar mejor el ángulo que forman los vectores del gradiente con el mallado pedido.
4. Campo de vectores en el instante inicial.
Dibujar el campo de vectores en los puntos del mallado del sólido, en [math]t=0[/math].
% Dibujar el campo de vectores en el mallado
x=0:0.2:10;
y=0:0.2:2;
%Realizamos el mallado:
[xx,yy]=meshgrid(x,y);
%Indicamos las componentes del vector desplazamiento:
Ux=2/5.*sin(xx);
Uy=Ux.*0;
%Dibujamos el campo de desplazamientos:
hold on
mesh(xx,yy,0*xx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
quiver(xx,yy,Ux,Uy,'r');
5. Placa antes y después de las deformaciones
Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec{u}[/math] (en [math]t=0[/math]). Dibujar ambos en la misma figura usando el comando [math]subplot[/math].
x=0:0.2:10;
y=0:0.2:2;
%Realizamos el mallado:
[xx,yy]=meshgrid(x,y);
%Indicamos las componentes del vector desplazamiento:
Ux=2/5.*sin(xx);
Uy=Ux.*0;
%Obtenemos las nuevas posiciones del mallado
%sumando el campo obtenido en el ejercicio 4
RUx=xx+Ux;
RUy=yy+Uy;
%Dibujamos la placa original
subplot(2,1,1)
mesh(xx,yy,0*xx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
%Dibujamos la placa tras haberla desplazado
subplot(2,1,2)
mesh(RUx,RUy,0*RUx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
Como se puede apreciar, al aplicar las tensiones y provocar las deformaciones, la placa se ensancha en la zona de la izquierda y a medida que nos desplazamos hacia la derecha, esta se encoge. Sin embargo, si nos seguimos desplazando hacia la derecha, vemos como el sólido se vuelve a ensanchar, para finalmente terminar encogiéndose de nuevo.
6. Divergencia del campo vectorial sobre la capa
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?
Conocido [math]\vec{u}[/math], calcularemos la divergencia:
[math]∇·\vec{u} = \frac{∂u1}{∂x}(2/5senx)+\frac{∂u2}{∂y}(0)+\frac{∂u3}{∂z}(0) = 2/5cosx[/math]
La divergencia de [math]\vec{u}[/math] es máxima en x=0 y en x=2π, es mínima en x=π y en x=3π y nula en x=π/2, en x=(3/2)*π y en x=(5/2)* π
%Dibujar las curvas de nivel de la divergencia
%Discretizamos las variables x e y
x=0:0.2:10;
y=0:0.2:2;
%Creamos el mallado
[xx,yy]=meshgrid(x,y);
%Calculamos la divergencia para cada punto
%del mallado:
DivU=(2/5)*cos(xx);
%Dibujamos la divergencia como un campo vectorial:
subplot(2,1,1)
hold on
mesh(xx,yy,0*xx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
quiver(xx,yy,DivU,DivU*0,'b');
%[M,c]=contour(xx,yy,DivU,[-0.4:0.05:0.4],'ShowText','on');
%c.LineWidth=2;
hold off
%Dibujamos las curvas del campo escalar definido por la divergencia de u:
subplot(2,1,2)
hold on
mesh(xx,yy,0*xx)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
[M,c]=contour(xx,yy,DivU,[-0.4:0.05:0.4],'ShowText','on');
c.LineWidth=2;
hold off
Como se comentó en el apartado anterior, y tal y como se puede apreciar en la figura, la divergencia nos da la razón con la deformación de la placa. Esto es porque cuando dicha divergencia es negativa, la placa se encoge, mientras que cuando la divergencia es positiva, esta se ensancha. Por lo que se puede ver, a la izquierda, las curvas de nivel son positivas (ensanchamiento) y a medida que nos movemos a la derecha van disminuyendo su valor hasta llegar a ser negativas (encogimiento), para luego volver a ser positivas y por último, volver a ser negativas.
7. Cálculo del rotacional
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?
Calcularemos el rotacional:
[math]∇ × \vec{u}= \begin{vmatrix}\vec{i} & \vec{j}& \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ u_{1}& u_{2}& u_{3} \end{vmatrix} = \begin{vmatrix} \vec{i} & \vec{j} & \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} &\frac{\partial }{\partial z} \\ 2/5senx & 0 & 0\end{vmatrix} = 0 [/math]
Por lo que el módulo será:
[math]|∇ × \vec{u}|= 0 [/math]
Al momento de calcular el rotacional, nos encontramos con que este es igual a 0. Esto, analíticamente hablando, quiere decir que el eje de rotación del sólido es nulo, por lo que este no gira sobre sí mismo y se mueve de manera uniforme.
8. Cálculo de tensiones normales en la dirección de los ejes principales
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).
% Dibujar las tensiones normales
x=0:0.2:10;
y=0:0.2:2;
% Realizamos el mallado:
[xx,yy]=meshgrid(x,y);
%Obtenemos las expresiones de las tensiones
TX=(6/5).*cos(xx);
TY=(2/5).*cos(xx);
TZ=(2/5).*cos(xx);
%Tensiones normales según el Eje X
subplot(3,1,1)
surf(xx,yy,TX)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
colorbar
caxis([0 0.75]);
title('Tensiones normales en la dirección del eje X')
% Tensiones normales según el Eje Y
subplot(3,1,2)
surf(xx,yy,TY)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
colorbar
caxis([0 0.75]);
title('Tensiones normales en la dirección del eje Y')
% Tensiones normales según el Eje Z
subplot(3,1,3)
surf(xx,yy,TZ)
axis equal
axis([-0.5 10.5 -0.5 2.5])
view(2)
colorbar
caxis([0 0.75]);
title('Tensiones normales en la dirección del eje Z')
9. Cálculo de 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.
Una vez realizados los cálculos, llegamos al resultado final de que las tensiones tangenciales respecto al plano ortogonal al vector i son nulas, por lo que, siguiendo instrucciones del enunciado, no consideramos necesaria la representación de ninguna tensión tangencial.
10. Tensiones de Von Mises
La tensión de Von Mises se define por la fórmula
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 el desarrollo del código, repetimos los procesos anteriores, creando ahora las variable MVonM y definiendo la función de Von Mises (VonMises). creamos un bucle en el que asignamos los valores de la tensión en cada punto y, por último, representamos la gráfica.
%Dibujar las tensiones de Von Mises
%Discretizamos las variables x e y
x=0:0.2:10;
y=0:0.2:2;
%Creamos el mallado
[xx,yy]=meshgrid(x,y);
%Creamos la matriz para recibir las tensiones de Von Mises en cada punto
%del mallado con la ayuda de dos bucles
[f c]=size(xx);
MVM=zeros(size(xx));
for i=1:f
for j=1:c
s1=(6/5)*cos(xx(i,j));
s2=(2/5)*cos(xx(i,j));
s3=(2/5)*cos(xx(i,j));
MVM(i,j)=sqrt(((s1-s2).^2+(s2-s3).^2+(s3-s1).^2)/2);
end
end
%Identificamos la máxima tensión de Von Mises
maxMVM=max(max(MVM));
%Identificamos las coordenadas de los puntos de máxima tensión de Von Mises
[YMVM,XMVM]=find(MVM==maxMVM);
hold on
%Dibujamos el mallado que representa los puntos de la placa
mesh(xx,yy,xx*0);
%Dibujamos las curvas de nivel de las tensiones de Von Mises
[M,c]=contour(xx,yy,MVM,'ShowText','on');
c.LineWidth=2;
axis equal
axis([-0.5 10.5 -0.5 2.5])
title('Tensiones de Von Mises')
%Dibujamos los puntos de máxima tensión
plot(x(XMVM),y(YMVM),'r*')
title('Representación de la tensión de Von Mises sobre la placa')
hold off
Al tratarse de un campo escalar, procederemos a representar los valores que adquiere dicha tensión tal y como se ha hecho con la temperatura, representando las curvas de nivel. Aparte de ello, se han indicado con asteriscos rojos los valores en los que la tensión de Von Mises es máxima. En este caso, los mayores valores de tensión se encuentran en la parte izquierda de la gráfica.
11. Cálculo de la velocidad
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
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´on 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} = 2/5\vec{i}[/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.
Calculamos [math]∇ · σ[/math]
Calculamos los campos vectoriales que se obtienen al hacer la divergencia de la matriz σ</math> dichos campos son:
[math]∇ · \vec{v1}[/math] = -6/5 sin(x-vt)
[math]∇ · \vec{v2}[/math] =0
[math]∇ · \vec{v3}[/math] =0
Sustituyendo en la ecuación de la elasticidad lineal y despejando la velocidad:
[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ= -2/5 v^2·sin(x-vt) \vec{i}+6/5 sin(x-vt) \vec{i} [/math]
despejando v
nos queda que V= 3^1/2
No tenemos ondas transversales.
12. Cálculo del módulo del desplazamiento vertical
Fijado ahora el punto [math](1/2, 1)[/math], calcular el módulo del desplazamiento vertical (dirección [math]\vec{j}[/math]) a lo largo del tiempo en el intervalo [math]t ∈ [0, 10][/math]
De acuerdo a los datos proporcionados por el enunciado, el vector u viene dado por la siguiente expresión:
Como se puede apreciar, dicho vector va en la dirección del eje X. Por lo que si calculamos los desplazamientos verticales pedidos en la dirección [math]\vec{j}[/math], nos saldría un vector nulo, y quedaría demostrado que no hay desplazamiento vertical en la dirección del eje Y.
