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

Trabajo realizado por estudiantes
Título	Grupo A5-Visualización de campos escalares y vectoriales en elasticidad.
Asignatura
Curso
Autores	María Rodríguez de Juan, Paula Sánchez García, Carlota Sancho Solanas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Tomaremos una placa plana que ocupa la mitad de un anillo circular centrado en el origen y comprendido entre los radios 1 y 2, en el plano y ≥ 0. Dibujar con los ejes en el cuadrado [−3, 3] × [−1, 3]. tomar como función de temperatura T(x, y) la que viene dada por:

[math] T(x,y) = x^2 + (y-1)^2 [/math],

1 Dibujo del sólido

Para comenzar dibujamos un mallado que represente los puntos interiores del sólido, tomando los ejes en el cuadrado [-3; 3] × [-1; 3] y como paso de muestreo h = 1/10 para las variables x e y.

%Mallado interior placa recatangular
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 mesh(xx,yy,yy*0);
 axis([-3,3,-1,3]);
view(2)

2 Dibujo curvas de nivel

Dibujamos también las curvas de nivel de la temperatura, que viene dada por el campo escalar : [math] T(x,y) = x^2 + (y-1)^2 [/math]. Según observamos en la gráfica, la temperatura varía de la zona azul (más fría) a la zona amarilla (más cálida), encontrando sus máximos en los puntos [-2,0] y [2,0].

u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 figure(1)
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 mesh(xx,yy,yy*0);
%Campo que define la temperatura
T=xx.^2+(yy-1).^2;
%Dibujar las curvas de nivel (contour)
subplot(2,1,1)
surf(xx,yy,T)
 axis([-3,3,-1,3]);
 view(2)
 title('Campo de temperaturas');
subplot(2,1,2)
contour(xx,yy,T)
 axis([-3,3,-1,3]);
 view(2)
 colorbar

3 Cálculo y dibujo del gradiente

Calculamos y dibujamos el gradiente observando en el dibujo que los vectores son ortogonales a las curvas de nivel ya que el gradiente muestra la dirección máxima de variación en cada punto.

[math]grad T = 2x\vec i\ + 2(y-1)\vec j\ [/math]

%Cálculo del gradiente
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 figure(1)
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 mesh(xx,yy,yy*0);
figure(1)
T=xx.^2+(yy-1).^2;
%Dibujar las curvas de nivel (contour)
contour(xx,yy,T)
view(2)
hold on
 colorbar
hold on
%representación del gradiente
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 figure(1)
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
mx=2*xx;
my=2*(yy-1);
quiver(xx,yy,mx,my);
axis equal 
axis([-3,3,-1,3]);
 view(2)
 title('Gradiente de T')

4 Dibujo campo de vectores

A continuación dibujamos el campo de vectores en los puntos del mallado del sólido, tomando el campo dado en el enunciado:

[math] \vec u(ρ,θ)=sen(θ)\frac{(ρ)}{10}(\vec e_{ρ}) + cos(θ)\frac{(ρ)}{5}(\vec e_{\theta}) [/math]

%Campo de vectores en el mallado del sólido
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 mesh(xx,yy,yy*0);
figure(1)
mx=yy/20;
my=yy/10;
mesh(xx,yy,yy*0);
hold on
quiver(xx,yy,mx,my);
 axis([-3,3,-1,3]);
view(2)
title('Campo U')

5 Dibujo antes y después del desplazamiento

Seguidamente representamos el sólido antes y después del desplazamiento dado por el campo de vectores ū. Como podemos observar el campo de desplazamientos solo varía en la dirección de los ejes [math] \vec e_ρ\ [/math] y [math] \vec e_θ\ [/math] .

%Sólido antes y después del desplazamiento
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
%Posición final r(x,y)=r(0)(x,y)+U(x,y)
rx=xx+(yy/20);
ry=yy+(yy/10);
%Situación inicial
figure
subplot(1,2,1)
surf(xx,yy,yy*0)
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
title('Antes de desplazarse')
 axis([-3,3,-1,3]);;
view(2);
%Después del desplazamiento
subplot(1,2,2)
surf(rx,ry,ry*0)
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
title('Después de desplazarse')
 axis([-3,3,-1,3]);
view(2);

6 Cálculo de la divergencia

A continuación calculamos la divergencia, que es una medida del cambio de volumen local debido al desplazamiento:

[math]\nabla·\vec u(ρ,θ) =\frac{1}{ρ} \frac{\partial}{\partial ρ}({sen(θ)}{\frac{(ρ^2)}{10}})+\frac{1}{ρ} \frac{\partial}{\partial θ}({cos(θ)}{\frac{(ρ)}{5}})= \frac{(sen(θ))}{5}-\frac{(sen(θ))}{5}=0 [/math].

La divergencia mide el cambio de volumen inducido por el campo, cuando el cuerpo está en un medio elástico si se desplaza y la divergencia es cero quiere decir que no ha cambiado su volumen y por lo tanto sigue siendo el mismo. En este caso como obtenemos una divergencia igual a cero, el movimiento de las moléculas no afecta a su volumen (densidad constante). Esto se suele dar en líquidos incompresibles.

7 Cálculo y dibujo del rotacional

[math]\nabla×\vec u(ρ,θ) =\frac{1}{ρ}\left|\begin{matrix} \vec e_ρ & \vec ρe_θ & \vec e_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \ sen(θ)\frac{ρ}{10} & \ cosθ\frac{ρ^2}{5} & \ 0 \end{matrix}\right|=\frac{3cosθ}{10}\vec e_z [/math].

[math]|\nabla×\vec u(ρ,θ)| =\frac{3cosθ}{10}[/math]

Como podemos observar los puntos que sufren un mayor rotacional son aquellos representados en color amarillo, estos tendrán una mayor capacidad de giro ya que el rotacional nos indica la capacidad que tienen los distintos puntos de nuestra placa para girar sobre otro punto determinado.

%Rotacional del campo vectorial
 u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 d=(cos(teta).*(3/10));
 surf(xx,yy,d);
 colorbar
 axis([-3,3,-1,3]);
 view(2)
 title('Módulo del rotacional')

8 Dibujo tensiones normales

A continuación, analizaremos las tensiones a las que se ve sometida la placa en un medio elástico lineal, isótropo y homogéneo donde los desplazamientos permiten escribir el tensor de la siguiente forma:

[math]σ(x,y)=λ∇.(\vec u)1 + 2με[/math]

Sabiendo que ε es la parte simétrica del ∇\(\vec u\) y conociendo λ=μ=1 (debido al material), calculamos las tensiones normales y las representamos en el dibujo.

8.1 Cálculo del tensor de deformaciones y de tensiones

• [math]\frac{\partial \vec u}{\partial ρ}=(\frac{1}{10})senθ \vec e_ ρ + (\frac{ρ}{10})senθ \frac{\partial \vec e_ ρ}{\partial ρ} +(\frac{1}{5})cosθ\vec e_θ + (\frac{ ρ}{5})cosθ \frac{\partial \vec e_θ }{\partial ρ}=[/math]

[math]=(\frac{1}{10})senθ \vec e_ ρ + (\frac{1}{5})cosθ \vec e_θ [/math].

• [math]\frac{\partial \vec u}{\partial θ}=(\frac{ρ}{10})cosθ \vec e_ ρ + (\frac{ρ}{10})senθ \frac{\partial \vec e_ ρ}{\partial θ} +(\frac{ρ}{5})senθ\vec e_θ + (\frac{ ρ}{5})cosθ \frac{\partial \vec e_θ }{\partial θ}=[/math]

[math]=(-\frac{ρ}{10})cosθ \vec e_ ρ (-\frac{ρ}{10})senθ \vec e_θ [/math].

• [math]\frac{\partial \vec u}{\partial z}=(\frac{ρ}{10})senθ \frac{\partial \vec e_ ρ}{\partial z} +(\frac{ρ}{5})cosθ \frac{\partial \vec e_θ }{\partial z}=[/math]

[math]=0 [/math].

Calculamos [math]\nabla{\vec u(ρ,θ,z)}[/math] y [math](\nabla{\vec u(ρ,θ,z)})^t[/math].

• [math]\nabla{\vec u(ρ,θ,z)}= \left(\begin{matrix} (\frac{1}{10})senθ & (-\frac{1}{10})cosθ & 0 \\ (\frac{1}{5})cosθ & (-\frac{1}{10})senθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].

• [math](\nabla{\vec u(ρ,θ,z)})^t = \left(\begin{matrix} (\frac{1}{10})senθ & (\frac{1}{5})cosθ & 0 \\ (-\frac{1}{10})cosθ & (-\frac{1}{10})senθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].

Sabemos que [math]λ[/math] y [math]μ[/math] toman valor uno, sabiendo que la divergencia es cero, la ecuación se reduce a dos veces el tensor de deformaciones que a su vez está multiplicado por un medio por tanto el dos que multiplica se anula con el que divide y el resultado del tensor de tensiones sería:

• [math]e(\vec σ(ρ,θ,z))= \left(\begin{matrix} (\frac{1}{5})senθ & (\frac{1}{10})cosθ & 0 \\ (\frac{1}{10})cosθ & (-\frac{1}{5})senθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].

%Apartado 9:Tensiones normales según e ro
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 sigma=(1/5).*sin(teta);
 surf(xx,yy,sigma);
axis([-3,3,-1,3]);
colorbar
view(2)
title('Tensiones normales en la dirección del eje e ro')

%Apartado 9:Tensiones normales según e teta
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 sigma=(-1/5).*sin(teta);
 surf(xx,yy,sigma);
axis([-3,3,-1,3]);
colorbar
view(2)
title('Tensiones normales en la dirección del eje e teta')

Pese a que los desplazamientos se realizan en plano, las tensiones no tienen por qué ser planas y se pueden representar en las tres direcciones del espacio [math] \vec e_ρ\ [/math], [math] \vec e_θ\ [/math] y [math] \vec e_z\ [/math] . En este caso podemos ver que la tensión es más elevada en la dirección [math] \vec e_ρ\ [/math] que en la dirección [math] \vec e_θ\ [/math] por tanto es mayor en la dirección radial que en la perpendicular a los radios.

9 Cálculo tensiones tangenciales

Calculamos las tensiones tangenciales, que vienen dadas por: σ.*\(\vec e_{\rho}\)+(\(\vec e_{\rho}\)*σ*\(\vec e_{\rho}\))\(\vec e_{\rho}\)

[math]|σ.*(\vec e_{ρ})+((\vec e_{ρ})*σ*(\vec e_{ρ}))(\vec e_{ρ})|= |\frac{(sen(θ))}{5}(\vec e_{ρ}) -\frac{(sen(θ))}{5}(\vec e_{ρ})+ \frac{(cos(θ)}{10}(\vec e_{\theta})|=\frac{(cos(θ)}{10} [/math].

%Tensiones tangenciales al plano ortogonal a e ro
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 sigma=(1/10).*cos(teta);
 surf(xx,yy,sigma);
axis([-3,3,-1,3]);
colorbar
view(2)

10 Dibujo de la tensión de Von Mises

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

[math] \sigma_V =\sqrt{\frac{(\sigma_1-\sigma_2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2}{2}} = \sqrt{3(\frac{sen(θ)^2}{25}+\frac{cos(θ)^2}{100})} [/math]

Habiendo sacado anteriormente los autovalores de [math] \sigma [/math]: [math] \sigma_1= 0 [/math], [math] \sigma_2=\sqrt{\frac{(sen(θ)^2}{25}+\frac{(cos(θ)^2}{100}}[/math], [math] \sigma_3=-\sqrt{\frac{(sen(θ)^2}{25}+\frac{(cos(θ)^2}{100}}[/math]

%Tension de Von Mises
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 x=(3/25).*(sin(teta).^2);
 y=(3/100).*(cos(teta).^2);
 sigmavm=sqrt(x+y);
 surf(xx,yy,sigmavm);
 colorbar
 axis([-3,3,-1,3]);
title('Tensión de Von Mises')

Como podemos observar en el gráfico la tensión de Von Mises alcanza valor máximo en el punto (-0,06;1,99;0,35)en coordenadas cartesianas.

11 Cálculo y dibujo del campo de fuerzas

Para calcular el campo de fuerzas F que actúa sobre la placa se aproxima usando la ecuación de la elasticidad lineal [math]\vec F=\nabla·\vec σ [/math]

[math]\nabla·\vec σ(ρ) =\frac{1}{ρ} \frac{\partial}{\partial ρ}({sen(θ)}{\frac{(ρ)}{5}})+\frac{1}{ρ} \frac{\partial}{\partial θ}({\frac{(cos(θ)}{5}})=\frac{sen(θ)}{10ρ} [/math]

[math]\nabla·\vec σ(θ) =\frac{1}{ρ} \frac{\partial}{\partial ρ}({cos(θ)}{\frac{(ρ)}{10}})+\frac{1}{ρ} \frac{\partial}{\partial θ}({\frac{(-sen(θ)}{5}})=\frac{-cos(θ)}{10ρ} [/math]

[math] \vec F=\frac{(-sen(θ))}{10ρ}(\vec e_{ρ}) + \frac{(cos(θ))}{10ρ}(\vec e_{\theta}) [/math]

%Campo F
u=1:0.1:2;
v=0:0.1:pi;
[ro,teta]=meshgrid(u,v);
 xx=ro.*cos(teta);
 yy=ro.*sin(teta);
 mx=(-1/10.*ro).*sin(teta);
 my=(1/10.*ro).*cos(teta);
 quiver(xx,yy,mx,my);
 axis([-3,3,-1,3]);
 view(2);
 title('Campo de fuerzas F')

Comparando este campo de fuerzas con el anterior podemos observar que en el campo u se producía principalmente deformación en la dirección [math] \vec e_ρ\ [/math] , y en este campo la principal variación en el desplazamiento sería en la dirección de [math] \vec e_θ\ [/math] .