Visualización de campos escalares y vectoriales en elasticidad. G19

Trabajo realizado por estudiantes
Título	Deformaciones de una placa plana. Grupo 6-A
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Mario Del Amo, Lucía Lázaro, Juan Carlos Martin, Claudia Xiang Martín
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Un placa plana está definida por dos sectores de circunferencias concéntricas de radios 1 y 2 centradas en el origen. Dichos sectores están limitados por el plano \(y≥\frac{|x|}{2}\).

Esta placa es objeto de estudio para su comportamiento frente a una función de temperatura:

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

Como se puede apreciar, la función viene dada en coordenadas cartesianas.

Se dispone también de un campo [math]\vec u [/math] que, a diferencia de la función de temperatura, viene dado en coordenadas cilíndricas:

[math]\vec u(ρ,θ)=\frac{1}{2}senθ \vec e_ρ [/math]

Este campo se empleará para obtener información acerca de la placa.

Contenido

1 . Mallado de la placa plana
2 . Curvas de nivel y gradiente de la temperatura
- 2.1 Curvas de nivel
  - 2.1.1 Temperatura máxima
- 2.2 Gradiente
3 . Ley de Fourier
4 . Campo de vectores en mallado sólido
5 . Sólido antes y después del desplazamiento
6 . Divergencia del campo
- 6.1 Divergencias máxima y mínima
7 . Rotacional del campo
- 7.1 Rotacional máximo
8 . Tensor de tensiones
9 . Tensiones tangenciales respecto al plano ortogonal
10 . Tensión de Von Mises
- 10.1 . Tensión máxima

1 . Mallado de la placa plana

Nuestra placa consiste en un cuarto de anillo centrado en el origen definido por la siguiente serie de condiciones:

• El anillo esta comprendido entre los radios 1 y 2

• Se encuentra en el plano \(y≥\frac{|x|}{2}\)

Debido a la forma de la placa será más conveniente utilizar coordenadas cilíndricas. De este modo, la última condición tendrá que ser pasada de cartesianas a cilíndricas, transformándose en [math]\rho sen\theta \geq \frac{\left | \rho cos\theta \right |}{2}[/math]. La variable ρ coincide con el radio de la plataforma, perteneciendo al intervalo [math]\left [ 1,2 \right ] [/math]. Sin embargo para obtener θ tendremos que despejarla de la expresión [math]\rho sen\theta \geq \frac{\left | \rho cos\theta \right |}{2}[/math]. Obtenemos [math]\tan \theta \geq \frac{1}{2}[/math] y [math]\tan \theta \lt -\frac{1}{2}[/math], por lo que θ tendrá que estar contenido en el intervalo [math]\left [ \arctan \frac{1}{2},180-\arctan \frac{1}{2} \right ][/math]

Mallado de la placa plana

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
mesh(X,Y,0*X);
axis([-3,3,-1,3]);
title('Mallado');
xlabel('Eje X');
ylabel('Eje Y');
view(2)

2 . Curvas de nivel y gradiente de la temperatura

2.1 Curvas de nivel

A continuación utilizaremos la función de temperatura definida en la introducción y dibujaremos sus curvas de nivel. Las curvas representan el efecto de la temperatura sobre la placa, permitiéndonos distinguir entre las zonas más cálidas (representadas con colores más claros como puede ser el verde o el amarillo) y las zonas más frías (representadas con colores más oscuros como el azul).

Curvas de nivel de la placa

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
mesh(X,Y,0*X);
Temp=cos(X.^2)+sin((Y-1).^2);
hold on
subplot(1,2,1);
axis([-3,3,-1,3]);
surf(X,Y,Temp);
view(2)
title('Grafica 1');
axis equal
xlabel('Eje X');
ylabel('Eje Y');
colorbar
subplot(1,2,2);
contour(X,Y,Temp,50);
axis([-3,3,-1,3]);
title('Placa 2D');
axis equal
xlabel('Eje X');
ylabel('Eje Y');
contour(X,Y,Temp,50);
colorbar
hold off

2.1.1 Temperatura máxima

Empleando el código mostrado abajo se ha obtenido la temperatura máxima que alcanza la placa. Esta temperatura es de 1.8315ºC en el punto (0,2) de la gráfica, como se puede apreciar en la figura superior derecha.

tempmax=max(max(Temp))

2.2 Gradiente

El gradiente representa la dirección en la cual la temperatura aumenta en mayor medida. Debido a esto es necesario que el gradiente sea perpendicular a las curvas de nivel, como se puede observar en la imagen del gradiente aumentado. Para realizar el gradiente de la temperatura utilizaremos la siguiente fórmula:

[math]\nabla T(x,y) =\frac{∂T}{∂x}\vec i + \frac{∂T}{∂y}\vec j [/math]

De esta forma nuestro gradiente de nuestra función [math]T(x,y)=cos(x^2)+sen((y-1)^2)[/math] tendrá como resultado:

[math]\nabla T(x,y) = \sin (x^{2})2x\vec i+\cos ((y-1)^{2})2(y-1)\vec j[/math]

Gradiente del campo atravesado por la placa

Zoom del gradiente

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
%Pasamos de coordenadas cilíndricas a cartesianas
X=R.*cos(T);
Y=R.*sin(T);
mesh(X,Y,0*X);
%Gradiente
dX=2.*X.*sin(X.^2);
dY=-2.*(Y-1).*cos((Y-1).^2);
hold on
quiver(X,Y,dX,dY);
hold off
view(2)
axis([-3,3,-1,3]);
title('Gradiente');
xlabel('Eje X');
ylabel('Eje Y');

3 . Ley de Fourier

En la Ley de Fourier, la energía calorífica [math]\vec Q [/math] viaja de acuerdo a la fórmula [math] \vec Q=-k\nabla T[/math], donde κ es la constante de conductividad térmica de la placa que supondremos κ = 1.

En la figura se muestra la variación de [math]\vec Q [/math] a lo largo de la placa y el programa Matlab/Octave empleado para su obtención.

Campo vectorial

h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
Temp=cos(X.^2)+sin((Y-1).^2);
Qx=-sin(X.^2).*2.*X;
Qy=-cos((Y-1).^2).*2.*(Y-1);
figure(1);
mesh(X,Y,0*X);
contour(X,Y,Temp);
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar;
hold on
quiver(X,Y,Qx,Qy);
axis equal
axis([-3,3,-1,3]);
view(2);
title('Energía calorífica. Ley de Fourier');
hold off

4 . Campo de vectores en mallado sólido

Un campo de vectores representa la distribución espacial de una magnitud vectorial, en este caso el campo es [math]\vec u(ρ,θ)=\frac{1}{2}senθ \vec e_ρ [/math].

El gráfico muestra el campo de vectores de fuerza para t=0. Sin embargo, el campo [math]\vec u [/math] no depende del tiempo.

También se muestra el programa Matlab/Octave para la obtención del gráfico mencionado.

Campo de vectores en la malla

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
mesh(X,Y,0.*X);
view(2)
axis([-3,3,-1,3]);
hold on
i=0.5.*sin(T).*cos(T);
j=0.5.*sin(T).*sin(T);
quiver(X,Y,i,j);
title('Campo de vectores');
hold off

5 . Sólido antes y después del desplazamiento

La placa sólida se enfrenta a dicho campo de vectores fuerza provocando un desplazamiento de la misma.

En la primera figura se muestra la situación inicial y final del desplazamiento, obtenida en Matlab/Octave.

Campo de vectores en la malla

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
subplot(1,2,1)
M=mesh(X,Y,0.*X);
view(2)
set(M,'EdgeColor','g');
axis([-3,3,-1,3]);
title('Placa inicial');
subplot(1,2,2)
u=0.5.*sin(T).*cos(T);
v=0.5.*sin(T).^2;
U=X+u;
V=Y+v;
N= mesh(U,V,0.*U);
view(2)
set(N,'EdgeColor','b');
axis([-3,3,-1,3]);
title('Placa final');

En esta figura se vuelven a mostrar la situaciones inicial y final superpuestas para apreciar la diferencia entre ambas.

Se puede percibir que la placa es mayor y cambia su posición en la situación final (con repecto a la situación inicial) siguiendo las líneas del campo vectorial que se mostró en el apartado 4.

Gráfica conjunta de la malla antes y después del desplazamiento

%Comparación entre sólidos
clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
u=0.5.*sin(T).*cos(T);
v=0.5.*sin(T).^2;
U=X+u;
V=Y+v;
figure
axis([-3,3,-1,3]);
N= mesh(U,V,0.*U);
set(N,'EdgeColor','b');
view(2)
title('Comparación entre sólido final e inicial');
xlabel('Eje X');
ylabel('Eje Y');
hold on
M= mesh(X,Y,0.*X);
set(M,'EdgeColor','g');
hold off

6 . Divergencia del campo

La divergencia mide la diferencia de flujo entrante y saliente de un campo vectorial sobre una superficie.

Por tanto, cuanto mayor sea la divergencia, mayor es la rapidez con la que se conduce el flujo a través de la superficie. Si la divergencia tiende a 0, su incompresibilidad es mayor.

La divergencia a partir de coordenadas cilíndricas, como ocurre con el campo dado, se calcula empleando:

[math] \nabla\cdot\vec u = \frac{1}{\rho}(\frac{\partial(ρu_\rho)}{\partial \rho} + \frac{\partial (u_\varphi)}{\partial \varphi} + \frac{\partial (ρu_z)}{\partial z}) [/math]

Concretamente con el campo [math]\vec u(ρ,θ)=\frac{1}{2}senθ\vec e_ρ [/math] dado el cálculo da lugar a que la divergencia sea:

[math] \nabla\cdot\vec u = \frac{1}{ρ}\frac{\partial(ρ\frac{1}{2}senθ)}{\partial \rho} =\frac{1}{2ρ}senθ\ [/math]

En la siguiente figura se observa el programa Matlab/Octave utilizado para obtener el gráfico de los valores de la divergencia en el cuarto de anillo en función del campo [math]\vec u [/math].

Divergencia del campo

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
div=(1./2.*R).*sin(T);
surf(X,Y,div);
view(2);
axis([-3,3,-1,3]);
title('Divergencia del campo');
xlabel('Eje X');
ylabel('Eje Y');
colorbar

6.1 Divergencias máxima y mínima

Empleando los códigos escritos en la parte inferior del apartado, se obtiene la máxima y mínima divergencia, las cuales son 0.9957 y 0.2236, respectivamente. La divergencia máxima se alcanza en el punto (0,2) de la gráfica superior y el mínimo en el extremo inferior derecho de la placa.

divmax=max(max(div))
divmin=min(min(div))

7 . Rotacional del campo

El rotacional consiste en obtener la tendencia de un campo a producir una rotación alrededor de un punto de la superficie. El resultado obtenido es un vector que indica la dirección de rotación de la superficie en dicho punto. Los cálculos a realizar, en este caso, en coordenadas cilíndricas son:

[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} \\ \vec e_ρ & ρ\vec e_θ & \vec e_z \end{matrix}\right|[/math].

Para el caso del campo [math]\vec u(ρ,θ)=\frac{1}{2}senθ\vec e_ρ [/math] ya definido, el rotacional es:

[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} \\ \frac{1}{2}senθ & 0 & 0 \end{matrix}\right|= \frac{-1}{2ρ}cosθ\vec e_z[/math]

Sin embargo, se pide el módulo del rotacional, que este caso es:

[math]|∇ × \vec{u}|= \frac{1}{2ρ}cosθ[/math]

En la siguiente figura se observa el programa Matlab/Octave utilizado para obtener el gráfico de los valores del rotacional en el cuarto de anillo en función del campo [math]\vec u [/math].

Rotacional del campo

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
mesh(X,Y,0.*X);
rot=cos(T)./(2.*R);
subplot(1,2,1);
surf(X,Y,rot);
view(2);
axis([-3,3,-1,3]);
colorbar
title('Rotacional');
xlabel('Eje X');
ylabel('Eje Y');
subplot(1,2,2);
surf(X,Y,rot);
axis equal;
title('Rotacional 3D');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');

7.1 Rotacional máximo

El rotacional máximo, obtenido mediante el código inferior, es 0.4472 que se alcanza en la zona en la que la divergencia alcanza el mínimo.

rotmax=max(max(rot))

8 . Tensor de tensiones

Aunque el desplazamiento depende de [math]\theta[/math], hay tensiones en las tres direcciones. Para su cálculo se emplean las siguientes fórmulas:

[math]\epsilon (\vec{u})=\frac{\triangledown \vec{u}+\triangledown \vec{u}^t}{2}[/math], que es la parte simétrica del gradiente del campo [math]\vec u[/math].

[math] \sigma =\lambda \triangledown \cdot \vec{u}\cdot Id + 2\mu \epsilon [/math], siendo el tensor de tensiones descrito por los desplazamientos en un medio elástico lineal, isótropo y homogéneo. Siendo [math] Id [/math] es el tensor identidad, [math]\lambda[/math] y [math]\mu[/math] son los coeficientes de Lamé y sus valores son 1.

Los resultados son:

[math]\epsilon (\vec{u}) = \left(\begin{matrix} 0 & \frac{1}{4}cosθ & 0 \\ \frac{1}{4}cosθ & \frac{1}{2}senθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].

[math] \sigma[/math] = \begin{pmatrix} \frac{1}{2\rho}sen\theta & \frac{1}{2}cos\theta & 0\\ \frac{1}{2}cos\theta & sen\theta(\frac{1}{2\rho}+1) & 0\\ 0 & 0 & \frac{1}{2\rho}sen\theta \end{pmatrix}

Las tensiones son ortogonales al plano, son tensiones normales a [math]\vec e[/math], que se calculan aplicando el tensor [math](\vec e\otimes\vec e)·\sigma[/math], o lo que es lo mismo [math]\vec e·\sigma·\vec e [/math].

Siendo [math]\vec e_\rho=\begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix}[/math], [math]\vec e_\theta=\begin{pmatrix} 0\\ 1\\ 0 \end{pmatrix}[/math] y [math]\vec e_z=\begin{pmatrix} 0\\ 0\\ 1 \end{pmatrix}[/math], las tensiones normales son:

Para [math]\vec e_\rho = \vec e_\rho·\sigma·\vec e_\rho = \frac{1}{2\rho}sen(\theta)[/math]

Para [math]\vec e_\theta= \vec e_\theta·\sigma·\vec e_\theta = \frac{1}{2}sen(\theta)+sen(\theta)[/math]

Para [math]\vec e_z = \vec e_z·\sigma·\vec e_z = \frac{1}{2\rho}sen(\theta )[/math]

Tensiones normales en ρ, θ y z

clc, clear;
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);

trho=(1./(2.*R)).*(sin(T));
subplot(3,1,1);
surf(X,Y,trho);
axis equal
title('Tensión normal en ρ');

ttheta=((1/2).*sin(T))+sin(T);
subplot(3,1,2);
surf(X,Y,ttheta);
axis equal;
title('Tensión normal en θ');

tz=(1./(2.*R)).*(sin(T));
subplot(3,1,3);
surf(X,Y,tz);
axis equal
title('Tensión normal en z');

En la figura de la derecha se muestran las tensiones normales a cada componente además del programa Matlab/Octave empleado para la obtención de dicha figura

9 . Tensiones tangenciales respecto al plano ortogonal

En este apartado nos centraremos en las tensiones tangenciales respecto al plano ortogonal. Teniendo σ, el cual ha sido calculado en el apartado anterior, se pueden determinar las tensiones tangenciales. Únicamente quedarán representadas en la gráfica aquellas las cuales no sean nulas. El plano que se usa de referencia para calcular las tensiones tangenciales en este caso es [math]\vec i [/math], es decir

[math]|σ·\vec i-(\vec i·σ·\vec i)\vec i| = \left |\begin{pmatrix} \frac{1}{2\rho }sin\theta & \frac{1}{2} cos\theta & 0\\ \frac{1}{2} cos\theta & sin\theta (\frac{1}{\rho }+1) &0 \\ 0 & 0 & \frac{1}{2\rho }sin\theta \end{pmatrix}\begin{pmatrix} 1\\ 0\\ 0\end{pmatrix} -\frac{1}{2\rho }\begin{pmatrix} 1\\ 0\\ 0\end{pmatrix}\right |=\begin{pmatrix} \frac{1}{2\rho }sin\theta \\ \frac{1}{2}cos\theta \\ 0 \end{pmatrix} -\begin{pmatrix} \frac{1}{2\rho }sin\theta \\ 0\\ 0\end{pmatrix}=\begin{pmatrix} 0\\ \frac{1}{2}cos\theta \\ 0\end{pmatrix}=\frac{1}{2}cos\theta [/math]

Tensiones tangenciales

%Realizamos el mallado
 r=1:0.2:2;
 t=linspace(pi/4,3*pi/4,8);
 [rr,tt]=meshgrid(r,t);
 xx=rr.*cos(tt);
 yy=rr.*sin(tt);
 %A partir de la tensión tangencial calculada dibujamos su gráfica
 tang = (1/2.*cos(tt));
 subplot(2,1,1)
 surf(xx,yy,tang)
 axis equal
 xlabel('x')
 ylabel('y')
 zlabel('z')
 view(2)
 colorbar
 title('Tensión tangencial a la dirección del eje eρ');
 subplot(2,1,2)
 surf(xx,yy,tang)
 axis equal
 xlabel('x')
 ylabel('y')
 zlabel('z')
 view(3)
 colorbar
 title('Tensión tangencial a la dirección del eje eρ');

10 . Tensión de Von Mises

Tensión de Von Mises

La tensión de Von Mises mide el comportamiento plástico y la ductilidad de un material. La calculamos mediante la fórmula:

[math]\sigma_{VM}=\sqrt{\frac{\left (\sigma_{1}- \sigma_{2} \right )^{2}+\left (\sigma_{2}- \sigma_{3} \right )^{2}+\left (\sigma_{3}- \sigma_{1} \right )^{2}}{2}}[/math]

siendo [math]\sigma_{1},\sigma_{2},\sigma_{3}[/math], los autovalores calculados a continuación.

Para el calculo de los autovalores empleamos Matlab/Octave. Debemos crear la matriz de la tensión, posteriormente sustituimos los datos y mediante el comando "eig" logramos obtener los autovalores de la matriz. De esta matriz 3x3 obtenemos los autovalores que se encuentran en la diagonal y los sustituimos en la fórmula de la tensión de Von Mises, consiguiendo generar la siguiente gráfica.

clc,clear
h=2/10;
r=1:h:2;
t=atan(1/2):h:pi-atan(1/2);
[R,T]=meshgrid(r,t);
X=R.*cos(T);
Y=R.*sin(T);
for i=1:length(t)
  for j=1:length(r)
  rho=R(i,j);
  theta=T(i,j);
M=[sin(theta)./(2.*rho),0.5.*cos(theta),0;0.5.*cos(theta),sin(theta).*(1./(2.*rho)+1),0;0,0,sin(theta)./2.*rho];
[v,u]=eig(M);
VonMises(i,j)=sqrt(((u(1,1)-u(2,2))^2+(u(2,2)-u(3,3))^2+(u(3,3)-u(1,1))^2)/2)
end
end

surf(X,Y,VonMises)
title('Tensión de Von Mises');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
axis equal

10.1 . Tensión máxima

La tensión máxima es una magnitud escalar que se utiliza como indicador de cuando reacciona un material al iniciar un comportamiento plástico. Este indicador nos dirá por que parte de la placa se producirá la rotura debido a que se trata de la zona de mayor estrés.