Elasticidad en Campos Escalares y Vectoriales (Coordenadas Cilíndricas)

Trabajo realizado por estudiantes
Título	Elasticidad en Campos Escalares y Vectoriales (Coordenadas Cilíndricas) - Grupo 4-C
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Juan Casteres Cortiñas Jorge Martínez Rodríguez-Malo Manuel Martín-Oar Faustmann Sandra Fuzhen Rodríguez Ibáñez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

En el siguiente trabajo se expondrán los cálculos relacionados con el estudio de diferentes campos escalares y vectoriales, así como su influencia en un sólido curvado. Para ello se hará uso de coordenadas cartesianas y cilíndricas.

Se considera una placa rectangular plana (en dimensión 2) que ocupa la región [math](x, y) ∈ [−3, 3] × [−1, 3][/math]. También se suponen definidas las cantidades físicas de: la temperatura [math]T(x,y)=sin(x\cdot 2+(y−3)\cdot2)[/math] y el campo de desplazamientos: [math] u (ρ, θ) = log (3 − ρ) 2 cos(2θ)\vec{e_ρ} [/math] (Proporcionadas en coordenadas cartesianas y cilíndricas respectivamente).
Con estos datos, se procede a responder a las siguientes cuestiones:

1 Dibujo Lámina

Para el dibujo del mallado se tomarán la región dada: una parte de un anillo. Como paso de muestreo, tomaremos [math]\frac{1}{10}[/math], contrariamente a los [math]\frac{2}{10}[/math] que propone el enunciado, para las variables x e y. Esto es debido a que un paso de muestreo menor produce una malla con más puntos, y por lo tanto, con más detalle. Es por ello, que la elección de [math]\frac{1}{10}[/math] proporciona un buen balance entre detalle y claridad para todos los problemas planteados a lo largo de este artículo.
Para llevar a cabo el dibujo de la lámina (sólido), se ha empleado el programa informático MATLAB, debido a sus buenas prestaciones y sus múltiples facilidades a la hora de graficar problemas matemáticos.
Al tratarse de una geometría circular, se emplearán coordenadas cilíndricas. Sin embargo, MATLAB únicamente trabaja en coordenadas cartesianas, por lo que todos los resultados se expresarán, de forma parametrizada, en coordenadas cartesianas.

Representación gráfica de una lámina sólida en forma de cuarto de círculo

%Dibujo Lámina (sólido)
%% Datos y Variables
h=1/10; %Coordenadas r y t
r=1:h:2;
t=pi/4:h:(3*pi)/4;
%% Cálculo de coordenadas
[U,V]=meshgrid(r,t); %Generación de la retícula
X=U.*cos(V); %parametrización
Y=U.*sin(V);
%% Dibujo de la malla
mesh(X,Y,0*X);
axis([-3,3,-1,3]); % definimos ejes
title('Lámina');
xlabel('X');
ylabel('Y');
grid on;
view(2);

2 Cálculo de las curvas de nivel del campo de temperaturas

La temperatura del sólido viene dada por el siguiente campo escalar:

[math]T(x,y)=sin(x\cdot 2+(y−3)\cdot2)[/math]

Haciendo uso del comando "contour", se obtiene el siguiente gráfico de las curvas de nivel:

Representación gráfica del mapa de calor de la función temperatura T

%% Cálculo función de Temperatura
T=sin(X.^2+(Y-3).^2);
%% Gráfico Mapa Térmico
figure(1);
surf(X,Y,T);
title('Temperatura');
view(2);
axis equal;
axis([-3,3,-1,3]);
xlabel('X');
ylabel('Y');

2.1 Representación curvas de nivel de T

Las curvas de nivel son líneas que unen los puntos que cuentan con la misma altitud. Esto nos muestra cómo varía la función temperatura de una forma muy visual. Para representarlas en MATLAB se utiliza el comando "contour".
En la imagen de la derecha de pueden apreciar unas zonas blancas en la gráfica, estas se deben a que ahí la temperatura es nula.

Representación gráfica de las curvas de nivel de la función temperatura T

%% Gráfico Curvas de Nivel
figure(2);
contour(X,Y,T,50);
title('Curvas de Nivel');
axis equal
axis([-3,3,-1,3]);
xlabel('X');
ylabel('Y');
colorbar

2.2 Punto máximo de T

Para la función [math]T(x,y)[/math] cabe recalcar que los valores máximos y mínimos se alcanzan en el intervalo [math]x ∈ [−1,1][/math] lo que nos da a entender que se trata de una función seno. El punto máximo de la función temperatura se alcanza en el punto (x,y).
Para ello, se utilizará el comando "".

2.3 Cálculo de [math]\nabla T[/math]

Para visualizar el campo vectorial gradiente, empleamos el comando "gradient" y lo representamos gráficamente.

2.4 Ortogonalidad de [math]\nabla T[/math] respecto a las curvas de nivel

Representando el campo vectorial [math]\nabla T[/math] sobre las curvas de nivel, se puede apreciar que los vectores del gradiente son ortogonales a las curvas de nivel.

%% Gradiente de la Temperatura
[FX,FY]=gradient(T,h); %Cálculo del gradiente discretizando h
%Derivadas parciales
%Fx=(2.*X)/(X.^2+1);
%Fy=((2.*Y)-8)/((Y-4).^2+1);
figure(4);
contour(X,Y,T,50);
title('\nablaT');
xlabel('X');
ylabel('Y');
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar
hold on    
quiver(X,Y,FX,FY)
axis equal
axis([-3,3,-1,3]);
view(2)

3 Ley de Fourier aplicada a la energía calorífica [math]\vec{Q}[/math]

La energía calorífica es un término utilizado en el campo de la termodinámica, para referirse a aquella energía asociada al calor LO ACABO EN INGLES Según la Ley de Fourier, la energía calorífica [math]\vec{Q}[/math] viaja de acuerdo a la fórmula [math]\vec{Q} = −\kappa\nabla T[/math], donde κ es la constante de conductividad térmica de la placa que supondremos κ = 1.
Sabiendo que el gradiente de la temperatura representa la dirección de más velocidad de cambio de la función. El cual viene dado por:

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

Por tanto, queda:

[math]\nabla T(x,y) =2\cdot cos[2\cdot (x+y-3)] \vec i+ 2\cdot cos[2\cdot (x+y-3) \vec j[/math]

Una propiedad característica del gradiente es la ortogonalidad a las curvas de nivel que hemos representado en el apartado anterior.

%% Ley de Fourier
figure(5);
title('Q=-k\nablaT');
xlabel('X');
ylabel('Y');
view(2);
hold on    
%Representación campo vectorial Fourier
quiver(X,Y,-FX,-FY) 
axis equal
axis([-3,3,-1,3])
view(2)

4 Dibujo de un campo vectorial sobre el mallado, en t=0

En primer lugar, tomamos el campo vectorial [math]\vec{u}(ρ, θ) = log (3 − ρ) 2 cos(2θ)\vec{e_ρ} [/math] y lo parametrizamos como se muestra en las líneas 3 y 5 del código. Esto lo hacemos para tener [math]\vec{u}[/math] expresado en coordenadas cartesianas, ya que MATLAB no puede representarlo expresado en cilíndricas.

Posteriormente, representamos dicho campo vectorial sobre nuestro sólido, obteniendo el gráfico que se muestra a continuación.

%%Campo en t=0
%Componente x campo vectorial (Nueva parametrización)
ui=((log(3-U))./2).*cos(2.*V).*cos(V);
%Componente y campo vectorial
uj=((log(3-U))./2).*cos(2.*V).*sin(V);
figure(6);
title('Campo de vectores u en t=0');
xlabel('X');
ylabel('Y');
view(2);
hold on
%Representación campo vectorial u
quiver(X,Y,ui,uj)
axis equal
axis([-3,3,-1,3])
view(2)

5 Desplazamiento provocado por el campo vectorial [math]\vec{u}[/math]

Dado el campo: [math] u (ρ, θ) = log (3 − ρ) 2 cos(2θ)\vec{e_ρ} [/math] queremos representar el desplazamiento que experimenta nuestro sólido a causa de él.

5.1 Dibujo antes del desplazamiento

Previo al desplazamiento, la posición del sólido es idéntica a la representada en el apartado 1.

%Inicio
subplot(2,2,1);
j=mesh(X,Y,X*0);
set(j,'EdgeColor','g')
title('Situación antes del desplazamiento');
xlabel('X');
ylabel('Y');
axis equal;
axis([-3,3,-1,3])
view(2);

5.2 Dibujo después del desplazamiento

Tras la actuación del campo vectorial [math]\vec{u}[/math], la posición del sólido es la siguiente:

%Después del desplazamiento
subplot(2,2,2);
mesh(X+ui,Y+uj,X*0);
title('Situación después del desplazamiento');
xlabel('X');
ylabel('Y');
axis equal;
axis([-3,3,-1,3])
view(2);

5.3 Comparación pre y post deplazamiento

%Comparación entre ambas situaciones
subplot(2,2,3.5);
mesh(X+ui,Y+uj,X*0);
hold on
j=mesh(X,Y,X*0);
set(j,'EdgeColor','g');
title('Comaparación');
xlabel('X');
ylabel('Y');
axis equal;
axis([-3,3,-1,3])
view(2);

6 Dibujo de [math]\nabla \cdot \vec{u}[/math] en t=0

El operador diferencial de la divergencia de campos vectoriales en cilíndricas viene dado por la fórmula de:

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

Tomando nuestro campo vectorial: [math] u (\rho, \theta) = log (3 − \rho) 2 cos(2\theta)\vec{e_\rho} [/math]

[math]\nabla\cdot\vec{u}=\frac{1}{ρ}\cdot[\frac{\partial}{\partial\rho}(\rho· log (3 − ρ) 2 cos(2θ))+\vec 0 ({\vec u_\theta })+\vec 0({\rho·\vec u_z })] = \frac{(-cos(2V)}{2 \cdot (3-U)}+ \frac{(log(3-U)\cdot cos(2V)}{2U} [/math]

%% Divergencia
figure(8)
div=((-cos(2.*V))./(2.*(3-U)))+((log(3-U).*cos(2.*V))./2.*U);
surf(X,Y,div)
axis equal 
axis([-3,3,-1,3]);
xlabel('X');
ylabel('Y');
title('divergencia')
colorbar
view(2)

7 Cálculo y dibujo de [math]|\nabla \times \vec{u}|[/math] en t=0

El operador diferencial del rotacional de campos vectoriales en cilíndricas viene dado por:

[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 u_ρ & \vec ρu_θ & \vec u_z \end{matrix}\right| = \frac{1}{ρ}\left|\begin{matrix} \vec e_ρ & \vec e_θ & \vec e_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \ log (3 − ρ) 2 cos(2θ) \vec{u_ρ} & 0 & 0 \end{matrix}\right| = (\frac {1}{U} \cdot (-sin(2\cdot V))\cdot log(3-U))[/math].

Donde los puntos con mayor rotacional quedan tal que:

figure(9)
ROT=abs((1./U).*(-sin(2*V)).*log(3.-U)); %Ez
subplot(1,2,1)
%Rotacional en 2D
pcolor(X,Y,ROT);
colorbar;
title('Rotacional en 2D');
xlabel('Eje X');
ylabel('Eje Y');
axis on;
axis([-3,3,-1,3]);
%Rotacional en 3D
subplot(1,2,2)
surf(X,Y,ROT);
view(3)
colorbar;
title('Rotacional en 3D');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
axis on
axis([-3,3,-1,3]);

8 Tensor Tensiones

En un medio elástico lineal, isótropo y homogéneo se puede describir el tensor de tensiones mediante la expresión:

[math]\sigma=\lambda\nabla\cdot\vec{u}1+2\mu\epsilon[/math]

Donde 1 es el tensor identidad, [math]\lambda=\mu=1[/math] y [math]\epsilon[/math] es el tensor deformaciones, que viene dado por [math]\epsilon(\vec{u})=(\nabla\vec{u}+\nabla\vec{u}^t)/2[/math]. Por tanto, si calculamos [math]\epsilon[/math] y [math]\sigma[/math] matricialmente, podemos obtener los siguientes apartados. La matriz [math]\epsilon[/math] se define como la parte simétrica del gradiente del campo [math]\vec{u}[/math], y su expresión matricial es la siguiente:
[math]\nabla\vec{u}=\begin{pmatrix} \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)} & \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{2\cdot\rho} & 0\\ \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{2\cdot\rho} & \frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & 0\\ 0 & 0 & 0 \end{pmatrix} [/math]
La matriz identidad multiplicada por la divergencia del campo [math]\vec{u}[/math] se expresa de la siguiente forma:
[math]\nabla\cdot\vec{u}1=\begin{pmatrix} \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & 0 & 0\\ 0 & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & 0\\ 0 & 0 & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} \end{pmatrix} [/math]
Por lo tanto:
[math]\sigma=\begin{pmatrix} \frac{3}{2}\cdot\frac{-cos(2\cdot\theta)}{3-\rho}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & 0\\ \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{3}{2}\cdot\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{\rho} & 0\\ 0 & 0 & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} \end{pmatrix} [/math]

%% Apartado 8
i_sig_i=(-3./2).*(cos(2.*V)./(3-U))+(log(3-U).*cos(2.*V))./(2.*U);
j_sig_j=(-cos(2.*V)./(2.*(3-U)))+(3./2).*((log(3-U).*cos(2.*V))./U);
k_sig_k=(-cos(2.*V)./(2.*(3-U)))+((log(3-U).*cos(2.*V))./(2.*U));
figure(10)
subplot(2,2,1)
surf(X,Y,i_sig_i)
view(2)
axis equal 
axis([-3,3,-1,3]);
colorbar
title('i·\sigma·i')
subplot(2,2,2)
surf(X,Y,j_sig_j)
view(2)
axis equal 
axis([-3,3,-1,3]);
colorbar
title('j·\sigma·j')
subplot(2,2,3.5)
surf(X,Y,k_sig_k)
view(2)
axis equal 
axis([-3,3,-1,3]);
colorbar
title('k·\sigma·k')

8.1 Tensiones normales en la dirección de [math]\vec{e_\rho}[/math]

Vienen definidas por [math]\vec{e_\rho}·\sigma·\vec{e_\rho}[/math]
[math]\vec{e_\rho}·\sigma·\vec{e_\rho}=\begin{pmatrix} 1 & 0 & 0 \end{pmatrix}\begin{pmatrix}\frac{3}{2}\cdot\frac{-cos(2\cdot\theta)}{3-\rho}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & 0\\ \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{3}{2}\cdot\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{\rho} & 0\\ 0 & 0 & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}[/math]

[math]\begin{pmatrix} \frac{3}{2}\cdot\frac{-cos(2\cdot\theta)}{3-\rho}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & 0 \end{pmatrix}\begin{pmatrix} 1 & 0 & 0 \end{pmatrix}=[/math][math]\frac{3}{2}\cdot\frac{-cos(2\cdot\theta)}{3-\rho}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho}[/math]

8.2 Tensiones normales en la dirección de [math]\vec{e_\theta}[/math]

Vienen definidas por [math]\vec{e_\theta}·\sigma·\vec{e_\theta}[/math]
[math]\vec{e_\theta}·\sigma·\vec{e_\theta}=\begin{pmatrix} 0 & 1 & 0 \end{pmatrix}\begin{pmatrix}\frac{3}{2}\cdot\frac{-cos(2\cdot\theta)}{3-\rho}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & 0\\ \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{3}{2}\cdot\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{\rho} & 0\\ 0 & 0 & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} \end{pmatrix}\begin{pmatrix} 0 & 1 & 0 \end{pmatrix}[/math]
[math]\begin{pmatrix} \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{3}{2}\cdot\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{\rho} & 0 \end{pmatrix}\begin{pmatrix} 0 & 1 & 0 \end{pmatrix}=[/math][math]\frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{3}{2}\cdot\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{\rho}[/math]

8.3 Tensiones normales en la dirección de [math]\vec{e_z}[/math]

Vienen definidas por [math]\\vec{e_z}·\sigma·\vec{e_z}[/math]
[math]\vec{e_z}·\sigma·\vec{e_z}=\begin{pmatrix} 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} \frac{3}{2}\cdot\frac{-cos(2\cdot\theta)}{3-\rho}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} & \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & 0\\ \frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho} & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{3}{2}\cdot\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{\rho} & 0\\ 0 & 0 & \frac{-cos(2\cdot\theta)}{2\cdot(3-\rho)}+\frac{log(3-\rho)\cdot cos(2\cdot\theta)}{2\cdot\rho} \end{pmatrix}\begin{pmatrix} 0 & 0 & 1 \end{pmatrix}[/math]
[math]\begin{pmatrix} 0 & 0 & 0 \end{pmatrix}\begin{pmatrix} 0 & 0 & 1 \end{pmatrix}=[/math][math]\frac{-log(3-\rho)\cdot sin(2\cdot\theta)}{\rho}[/math]

9 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math]

%% Apartado 9
tti=sqrt(((log(3-U).*sin(2.*V).*sin(V)./U).^2))
ttj=sqrt(((-log(3-U).*sin(2.*V).*cos(V)./U).^2));
figure(11)
hold on
mesh(X,Y,X.*0)
quiver(X,Y,tti,ttj)
axis equal 
axis([-3,3,-1,3]);
title('tensiones tangenciales respecto al plano ortogonal a i');

10 Tensión de Von Mises

La tensión de Von Mises viene definida por la fórmula:

[math] \sigma_{VM}= \sqrt{\frac{(\sigma_1-\sigma_1)^2+(\sigma_2-\sigma_3)^2(\sigma_3-\sigma_1)^2}{2}}[/math]

donde [math]\sigma_1 , \sigma_2 , \sigma_3 [/math] son los autovalores de σ (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).
Para sacar los autovalores en MATLAB se utiliza el comando eig.m