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≥|x|\).

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.

1 . Mallado de la placa plana

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

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

Gradiente del campo atravesado por 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);
%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 Q viaja de acuerdo a la fórmula Q = −κ∇T, donde κ es la constante de conductividad térmica de la placa que supondremos κ = 1.

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

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

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('Grafica 1');
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('Grafica 2');

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

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');

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 quedaran representadas en la gráfica aquellas las cuales no sean nulas. El plano que se usa de referencia en este caso es [math]\vec i [/math], es decir [math]|σ·\vec i-(\vec i·σ·\vec i)\vec i| [/math]