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 una placa plana que ocupa un cuarto de anillo circular centrado en el origen y comprendido entre los radios 1 y 2. Dicha placa está contenida en el plano:
Este plano limita la placa circular entre [math]\frac{\pi}{8}[/math] y [math]\frac{7\pi}{8}[/math] como podemos observar en la siguiente figura.
En todas las figuras que se mostrarán a lo largo del artículo, dibujaremos los ejes en la región [math](x, y) ∈ [−3, 3] × [−1, 3][/math].
También se suponen definidas las siguientes cantidades físicas:
- 1. El campo escalar temperatura: [math]T(x,y)=sin(x^2+(y−3)^2)[/math], expresado en coordenadas cartesianas
2. El campo vectorial de desplazamientos: [math] \vec{u}(ρ, θ) =\frac {log (3 − ρ)} {2}\cdot cos(2θ)\vec{e_ρ} [/math], expresado en coordenadas cilíndricas
Con estos datos, se procede a responder a las siguientes cuestiones:
Contenido
- 1 Dibujo Lámina
- 2 Cálculo de las curvas de nivel del campo de temperaturas
- 3 Ley de Fourier aplicada a la energía calorífica [math]\vec{Q}[/math]
- 4 Dibujo de un campo vectorial sobre el mallado, en t=0
- 5 Desplazamiento provocado por el campo vectorial [math]\vec{u}[/math]
- 6 Dibujo de [math]\nabla \cdot \vec{u}[/math] en t=0
- 7 Cálculo y dibujo de [math]|\nabla \times \vec{u}|[/math] en t=0
- 8 Tensor Tensiones
- 9 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{e_\rho}[/math]
- 10 Tensión de Von Mises
1 Dibujo Lámina
Para el dibujo del mallado se tomará la región dada: una parte de un anillo. Como paso de muestreo para las variables [math]x[/math] e [math]y[/math], tomaremos [math]\frac{1}{10}[/math], contrariamente a los [math]\frac{2}{10}[/math] que propone el enunciado. 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 parametrizarán para expresarlos de tal forma.
%Dibujo Lámina (sólido)
%% Datos y Variables
h=1/10; %Coordenadas r y t
r=1:h:2;
t=pi/8:h:(7*pi)/8;
%% 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:
Si introducimos la expresión del campo en MATLAB, podemos representarlo gráficamente con el comando "surf", obteniendo la imagen que se presenta a continuación.
%% 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 del mismo valor. Esto nos muestra de una forma muy visual de la variación de la función temperatura. Es posible representar estas curvas en MATLAB mediante el comando "contour".
Las zonas de la gráfica que presentan colores más azulados son aquellas en las que la temperatura es menor, mientras que en las más amarillentas la temperatura es mayor.
Las partes en las que aparecen zonas blancas son aquellas en las que la temperatura es nula.
%% 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 trigonométrica seno.
En los dos gráficos obtenidos, tanto el del campo escalar temperatura como el de las curvas de nivel de dicho campo, podemos hacer zoom y observar que en ambos casos que el punto de la gráfica en el que la temperatura es mayor es aproximadamente el (0,2).
2.3 Cálculo de [math]\vec{\nabla T}[/math]
Podríamos calcular el gradiente del campo escalar T(x,y) analíticamente mediante las fórmulas vistas en clase. Sin embargo, introduciendo el campo escalar "T" y el paso "h" utilizado en apartados anteriores, el comando "gradient" nos lo calcula directamente.
El comando "contour" nos permite dibujar las curvas de nivel (apartado 2.1) y el comando "quiver" nos permite representar el campo vectorial gradiente. Por tanto, si representamos ambos a la vez mediante el comando "hold on", obtenemos el gráfico de [math]\vec{\nabla T}[/math] sobre las curvas de nivel de [math]T(x,y)[/math].
En esta imagen podemos apreciar que los vectores del gradiente son perpendiculares a las curvas de nivel. Esto se debe a que el gradiente es la dirección de máxima variación de la temperatura, mientras que en las curvas de nivel no hay variación de temperatura. Por tanto, el gradiente debe ser perpendicular a las curvas de nivel.
%% Gradiente de la Temperatura
[FX,FY]=gradient(T,h); %Cálculo del gradiente discretizando h
%Derivada parcial respecto de X=Fx=2.*cos(2.*X+2.*(Y-3));
%Derivada parcial respecto de Y=Fy=2.*cos(2.*X+2.*(Y-3));
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 que recibe, desprende o transmite un cuerpo. Un ejemplo físico muy común de la aplicación de esta fórmula es la transmisión de calor, y por tanto, variación de temperatura a lo largo de un sólido.
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 [math]\kappa[/math] 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:
El motivo por el cual la ley de Fourier incorpora un signo menos, tiene una explicación meramente física. Por el 2º Principio de la Termodinámica sabemos que al juntar un foco caliente y uno frío, el calor se transmite del foco caliente al foco frío.
Por tanto, queda:
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
Partimos del campo vectorial dado:
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 si está expresado en cilíndricas.
Posteriormente, representamos dicho campo vectorial sobre nuestro sólido mediante la función "quiver", 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]\vec{u}(ρ, θ) = \frac {log (3 − ρ)} {2} \cdot cos(2θ)\vec{e_ρ} [/math] queremos representar el desplazamiento que experimenta nuestro sólido a causa de él.
- 1. Imagen de arriba a la izquierda: situación del sólido antes del desplazamiento. Coincide con la situación representada en el apartado 1.
- 2. Imagen de arriba a la derecha: situación del sólido después del desplazamiento. Observamos cómo ciertos puntos de la placa se han desplazado hacia abajo, justo la dirección a la que apuntan los vectores del campo vectorial representado en el apartado 4.
- 3. Imagen de abajo: imágenes 1 y 2 superpuestas. Se puede apreciar mejor la diferencia entre las dos situaciones.
%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);
%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);
%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
La divergencia mide la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial sobre la superficie que rodea a un volumen de control.
En coordenadas cilíndricas, el operador diferencial de la divergencia de campos vectoriales viene dado por:
Tomando nuestro campo vectorial: [math]\vec{u}(ρ, θ) = \frac {log (3 − ρ)} {2} \cdot cos(2θ)\vec{e_ρ} [/math], sustituimos en la ecuación de arriba:
Siendo U y V los parámetros asociados respectivamente a las coordenadas [math]\rho[/math] y [math]\theta[/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)6.1 Análisis de [math]\nabla \cdot \vec{u}[/math]
Se pide determinar de forma analítica los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. Por lo que habiendo calculado el valor de [math]\nabla\cdot\vec{u}[/math] en el apartado anterior, se procede a utilizar la matriz Hessiana y sacar las derivadas del campo:
6.2 Representación del cambio de volumen
7 Cálculo y dibujo de [math]|\nabla \times \vec{u}|[/math] en t=0
El rotacional de un campo vectorial muestra la tendencia de un campo a inducir rotación alrededor de un punto. En coordenadas cilíndricas, el operador diferencial del rotacional de campos vectoriales viene dado por:
Sustituyendo para los valores de nuestro campo y operando, obtenemos:
Tras calcular [math]\nabla \times \vec{u}[/math], sabemos que es un campo vectorial. Sin embargo, lo pedido es [math]|\nabla \times \vec{u}|[/math], que es un campo escalar. Por tanto, se hace se procede a hacer el módulo del rotacional calculado anteriormente.
Al introducir el campo en MATLAB, parametrizamos poniéndolo en función de U y V:
Por tanto, el módulo será:
Si representamos gráficamente el campo escalar calculado, obtenemos las figuras que se muestran debajo.
figure(9)
ROT=sqrt(((1./U).*(-sin(2*V)).*log(3.-U)).^2); %Módulo del campo escalar.
%%ROT=abs((1./U).*(-sin(2*V)).*log(3.-U)); %Valor absoluto del campo.
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]);
hold on
plot(0.7,0.7,'ro') %Máximo 1 en gráfico 2D
hold on
plot(-0.65,0.75,'ro') %Máximo 2 en gráfico 2D
%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]);
hold on
plot3(0.7,0.7,0.693,'ro') %Máximo 1 en gráfico 3D
hold on
plot3(-0.65,0.75,0.693,'ro') %Máximo 2 en gráfico 3D
7.1 Cálculo punto máximo del rotacional
Como puntos máximos de rotacional, como se puede observar en el gráfico, nos quedan los siguientes puntos:
Los cuales están señalados en ambos gráficos con un círculo rojo.
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:
Donde 1 es el tensor identidad, [math]\lambda=\mu=1[/math] , [math]\epsilon[/math] es el tensor deformaciones que viene dado por [math]\epsilon(\vec{u})=(\nabla\vec{u}+\nabla\vec{u}^t)/2[/math] y [math]\nabla\cdot\vec{u}[/math] es la divergencia del campo [math]\vec{u}[/math].
Partiendo de la expresión matricial del gradiente, podemos calcular su parte simétrica, es decir, [math]\epsilon[/math]:
La matriz identidad multiplicada por la divergencia del campo [math]\vec{u}[/math] da como resultado la matriz:
Por lo tanto, operando con las matrices anteriores se obtiene la matriz del tensor tensiones:
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]. Por tanto, multiplicando la matriz [math]\sigma[/math] a ambos lados por el vector unitario [math]\vec{i}[/math], obtenemos las tensiones normales en la dirección de [math]\vec{i}[/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]. Por tanto, multiplicando la matriz [math]\sigma[/math] a ambos lados por el vector unitario [math]\vec{j}[/math], obtenemos las tensiones normales en la dirección de [math]\vec{j}[/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]. Por tanto, multiplicando la matriz [math]\sigma[/math] a ambos lados por el vector unitario [math]\vec{k}[/math], obtenemos las tensiones normales en la dirección de [math]\vec{k}[/math].
8.4 Código Matlab
%% 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')
Los tres campos escalares resultantes representan las tensiones normales en la dirección que marca el eje [math]\vec{e_\rho}[/math], las tensiones normales en la dirección que marca el eje [math]\vec{e_\theta}[/math], y las correspondientes al eje [math]\vec{e_z}[/math].
9 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{e_\rho}[/math]
Dichas tensiones vienen definidas por la expresión:
Puesto que [math](\vec{e_\rho}\cdot\sigma\cdot\vec{e_\rho})[/math] lo tenemos calculado del apartado 8.2, bastará con multiplicar el valor hallado por [math]\vec{e_\rho}[/math].
[math]\sigma\cdot\vec{e_\rho}[/math] se calcula de la siguiente forma:
Por lo tanto, el vector solución resultante de la resta de los anteriores es:
Dado que en Matlab se trabaja en coordenadas cartesianas, Es necesario pasar de la base física cilíndrica a la cartesiana.
Por lo tanto el módulo queda tal que:
tti=((log(3-U).*sin(2.*V).*sin(V)./U));
ttj=((-log(3-U).*sin(2.*V).*cos(V)./U));
mod=sqrt((tti.^2)+(ttj.^2));
figure(11)
hold on
mesh(X,Y,X.*0)
surf(X,Y,mod)
axis equal
axis([-3,3,-1,3]);
title('módulo de 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:
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
%% Matriz Tensión de Von Misses
sig=[];
M_VonMises=zeros(length(t),length(r));
a=@(UU,VV)(-3./2).*(cos(2.*VV)./(3-UU))+(log(3-UU).*cos(2.*VV))./(2.*UU);
b=@(UU,VV)(-cos(2.*VV)./(2.*(3-UU)))+(3./2).*((log(3-UU).*cos(2.*VV))./UU);
c=@(UU,VV)(-cos(2.*VV)./(2.*(3-UU)))+((log(3-UU).*cos(2.*VV))./(2.*UU));
d=@(UU,VV)((-log(3-UU).*sin(2.*VV))./UU);
%fila1=[a,d,0];
%fila2=[d,b,0];
%fila3=[0,0,c];
%sigma=[fila1;fila2;fila3]
% Se ha de asignar a la matriz "M_VonMises" los valores
% de la tensión de Von Mises en cada punto de la placa.
for i=1:length(t)
for j=1:length(r)
sig(1,1)=a(U(i,j),V(i,j));
sig(1,2)=d(U(i,j),V(i,j));
sig(1,3)=0;
sig(2,1)=d(U(i,j),V(i,j));
sig(2,2)=b(U(i,j),V(i,j));
sig(2,3)=0;
sig(3,1)=0;
sig(3,2)=0;
sig(3,3)=c(U(i,j),V(i,j));
[v,u]=eig(sig);%Autovalores
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
% Von Mises.10.1 Representación de la Tensión de Von Mises
%Representación Von Mises
figure(12)
surf(X,Y,M_VonMises);
% Se pone título a la gráfica.
title('Von Mises');
% Se da nombre a los ejes.
xlabel('Eje de las X');
ylabel('Eje de las Y');
axis equal
axis([-3,3,-1,3]);
% Se da equidistancia a los ejes.
% Se le aplica una barra con colores en función
% de la tension de Von Mises en cada punto.
colorbar;
view(2)
figure(13)
surf(X,Y,M_VonMises);
% Se pone título a la gráfica.
title('Von Mises');
% Se da nombre a los ejes.
xlabel('Eje de las X');
ylabel('Eje de las Y');
axis equal
axis([-3,3,-1,3]);
% Se da equidistancia a los ejes.
% Se le aplica una barra con colores en función
% de la tension de Von Mises en cada punto.
colorbar;
view(3)
